ameriflux/11.ameriflux.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Thu, 26 Mar 2009 14:02:21 -0400
changeset 1 4be95183fbcd
permissions -rw-r--r--
Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
forrest@0
     1
;************************************************************
forrest@0
     2
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0
     3
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0
     4
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0
     5
;************************************************************
forrest@0
     6
procedure set_line(lines:string,nline:integer,newlines:string) 
forrest@0
     7
begin
forrest@0
     8
; add line to ascci/html file
forrest@0
     9
    
forrest@0
    10
  nnewlines = dimsizes(newlines)
forrest@0
    11
  if(nline+nnewlines-1.ge.dimsizes(lines))
forrest@0
    12
    print("set_line: bad index, not setting anything.") 
forrest@0
    13
    return
forrest@0
    14
  end if 
forrest@0
    15
  lines(nline:nline+nnewlines-1) = newlines
forrest@0
    16
;  print ("lines = " + lines(nline:nline+nnewlines-1))
forrest@0
    17
  nline = nline + nnewlines
forrest@0
    18
  return 
forrest@0
    19
end
forrest@0
    20
;*************************************************************
forrest@0
    21
begin
forrest@0
    22
forrest@0
    23
  plot_type     = "ps"
forrest@0
    24
  plot_type_new = "png"
forrest@0
    25
forrest@0
    26
;------------------------------------------------------
forrest@0
    27
; edit table.html of current model for movel1_vs_model2
forrest@0
    28
forrest@0
    29
 if (isvar("compare")) then
forrest@0
    30
    html_name2 = compare+"/table.html"  
forrest@0
    31
    html_new2  = html_name2 +".new"
forrest@0
    32
 end if
forrest@0
    33
forrest@0
    34
;------------------------------------------------------
forrest@0
    35
; edit table.html for current model
forrest@0
    36
forrest@0
    37
 html_name = model_name+"/table.html"  
forrest@0
    38
 html_new  = html_name +".new"
forrest@0
    39
forrest@0
    40
;------------------------------------------------------
forrest@0
    41
forrest@0
    42
  nmonth = 12
forrest@0
    43
forrest@0
    44
; for nee, gpp, and ar
forrest@0
    45
; observed unit is gC/m2/day
forrest@0
    46
; model    unit is gC/m2/s
forrest@0
    47
; to change to observed unit,
forrest@0
    48
forrest@0
    49
  factor_flux = 86400.
forrest@0
    50
  
forrest@0
    51
; for incident solar radiation,
forrest@0
    52
; observed Rg_f unit is MJ/m2/day
forrest@0
    53
; model (FSDS)  unit is  W/m2
forrest@0
    54
; to change to model unit,
forrest@0
    55
forrest@0
    56
  factor_rad = 1.e6/86400.
forrest@0
    57
forrest@0
    58
;************************************************
forrest@0
    59
; observed data info
forrest@0
    60
;************************************************
forrest@0
    61
forrest@0
    62
 station = (/"ARM_Oklahoma" \
forrest@0
    63
            ,"ARM_Oklahoma_burn" \
forrest@0
    64
            ,"ARM_Oklahoma_control" \
forrest@0
    65
            ,"Atqasuk" \
forrest@0
    66
            ,"Audubon" \
forrest@0
    67
            ,"AustinCary" \
forrest@0
    68
            ,"Bartlett" \
forrest@0
    69
            ,"Bondville" \
forrest@0
    70
            ,"Brookings" \
forrest@0
    71
            ,"Donaldson" \
forrest@0
    72
            ,"Duke_Forest_Hardwoods" \
forrest@0
    73
            ,"Duke_Forest_Open_Field" \
forrest@0
    74
            ,"Duke_Forest_Pine" \
forrest@0
    75
            ,"Fermi_Ag" \
forrest@0
    76
            ,"Fermi_Prairie" \
forrest@0
    77
            ,"Flagstaff_Managed" \
forrest@0
    78
            ,"Flagstaff_Unmanaged" \
forrest@0
    79
            ,"Flagstaff_Wildfire" \
forrest@0
    80
            ,"FortPeck" \
forrest@0
    81
            ,"FreemanRanch_mesquite" \
forrest@0
    82
            ,"Goodwin_Creek" \
forrest@0
    83
            ,"HarvardForest" \
forrest@0
    84
            ,"HarvardForestHemlock" \
forrest@0
    85
            ,"HowlandForestMain" \
forrest@0
    86
            ,"HowlandForestWest" \
forrest@0
    87
            ,"Ivotuk" \
forrest@0
    88
            ,"KendallGrasslands" \
forrest@0
    89
            ,"KennedySpaceCenterPine" \
forrest@0
    90
            ,"KennedySpaceCenterScrub" \
forrest@0
    91
            ,"LittleProspect" \
forrest@0
    92
            ,"LostCreek" \
forrest@0
    93
            ,"Mead-irrigated" \
forrest@0
    94
            ,"Mead-irrigated-rotation" \
forrest@0
    95
            ,"Mead-rainfed" \
forrest@0
    96
            ,"Metolius_2nd_YoungPonderosaPine" \
forrest@0
    97
            ,"MetoliusEyerly" \
forrest@0
    98
            ,"MetoliusIntermediatePine" \
forrest@0
    99
            ,"MetoliusOldPonderosaPine" \
forrest@0
   100
            ,"MissouriOzark" \
forrest@0
   101
            ,"Mize" \
forrest@0
   102
            ,"MorganMonroe" \
forrest@0
   103
            ,"NiwotRidge" \
forrest@0
   104
            ,"NorthCarolina_cc" \
forrest@0
   105
            ,"NorthCarolina_lp" \
forrest@0
   106
            ,"ParkFalls" \
forrest@0
   107
            ,"Rayonier" \
forrest@0
   108
            ,"SantaRita" \
forrest@0
   109
            ,"SkyOaks_Old" \
forrest@0
   110
            ,"SkyOaks_PostFire" \
forrest@0
   111
            ,"SkyOaks_Young" \
forrest@0
   112
            ,"SylvaniaWilderness" \
forrest@0
   113
            ,"Toledo" \
forrest@0
   114
            ,"Tonzi" \
forrest@0
   115
            ,"UCI_1850" \
forrest@0
   116
            ,"UCI_1930" \
forrest@0
   117
            ,"UCI_1964" \
forrest@0
   118
            ,"UCI_1964wet" \
forrest@0
   119
            ,"UCI_1981" \
forrest@0
   120
            ,"UCI_1989" \
forrest@0
   121
            ,"UCI_1998" \
forrest@0
   122
            ,"UMBS" \
forrest@0
   123
            ,"Vaira" \
forrest@0
   124
            ,"WalkerBranch" \
forrest@0
   125
            ,"WillowCreek" \
forrest@0
   126
            ,"WindRiver" \
forrest@0
   127
            ,"Wisconsin_ihw" \
forrest@0
   128
            ,"Wisconsin_irp" \
forrest@0
   129
            ,"Wisconsin_mrp" \
forrest@0
   130
            ,"Wisconsin_myjp" \
forrest@0
   131
            ,"Wisconsin_pb" \
forrest@0
   132
            ,"Wisconsin_rpcc" \
forrest@0
   133
            ,"Wisconsin_yhw" \
forrest@0
   134
            ,"Wisconsin_yjp" \
forrest@0
   135
            ,"Wisconsin_yrp" \
forrest@0
   136
            /)
forrest@0
   137
forrest@0
   138
 field   = (/"NEE Flux" \
forrest@0
   139
            ,"Shortwave Incoming" \
forrest@0
   140
            ,"Latent Heat" \
forrest@0
   141
            ,"Sensible Heat" \
forrest@0
   142
            ,"GPP Flux" \
forrest@0
   143
            ,"Respiration" \
forrest@0
   144
            /)
forrest@0
   145
forrest@0
   146
 field_unit = (/"gC/m2/day" \
forrest@0
   147
               ,"W/m2" \
forrest@0
   148
               ,"W/m2" \
forrest@0
   149
               ,"W/m2" \
forrest@0
   150
               ,"gC/m2/day" \
forrest@0
   151
               ,"gC/m2/day" \
forrest@0
   152
               /)
forrest@0
   153
forrest@0
   154
 nstation = dimsizes(station)
forrest@0
   155
 nfield   = dimsizes(field)
forrest@0
   156
forrest@0
   157
;========================================================================
forrest@0
   158
; get observed info: number of year, first/last year, lat, lon
forrest@0
   159
; and annual data
forrest@0
   160
forrest@0
   161
 dir_root = diro + "ameriflux/"
forrest@0
   162
forrest@0
   163
 year_station = new ((/nstation/),integer) ; number of year
forrest@0
   164
 year_ob      = new ((/nstation/),string)  ; observed year
forrest@0
   165
 year_ob_i    = new ((/nstation/),integer) ; first year
forrest@0
   166
 year_ob_f    = new ((/nstation/),integer) ; last year
forrest@0
   167
 lat_ob       = new ((/nstation/),float)   ; latitude
forrest@0
   168
 lon_ob       = new ((/nstation/),float)   ; longitude 
forrest@0
   169
forrest@0
   170
 data_ob_ann  = new ((/nfield, nmonth, nstation/),float)
forrest@0
   171
forrest@0
   172
 do n = 0, nstation-1
forrest@0
   173
forrest@0
   174
    dir_f = dir_root + station(n)+"/"
forrest@0
   175
    fil_f = "timeseries_L4_m.nc"
forrest@0
   176
    fo   = addfile (dir_f+fil_f,"r")
forrest@0
   177
 
forrest@0
   178
    lat_ob(n) = fo->lat
forrest@0
   179
    lon_ob(n) = fo->lon
forrest@0
   180
forrest@0
   181
    year      = fo->year
forrest@0
   182
forrest@0
   183
    year_station(n) = dimsizes(year)
forrest@0
   184
    year_ob_i(n)    = year(0)
forrest@0
   185
    year_ob_f(n)    = year(year_station(n)-1)
forrest@0
   186
    year_ob(n)      = year_ob_i(n) + "-" + year_ob_f(n)
forrest@0
   187
 
forrest@0
   188
    delete (year)
forrest@0
   189
 
forrest@0
   190
    data = fo->NEE_or_fMDS
forrest@0
   191
    data_ob_ann(0,:,n) = dim_avg(data(month|:,year|:))
forrest@0
   192
forrest@0
   193
    data = fo->Rg_f
forrest@0
   194
    data_ob_ann(1,:,n) = dim_avg(data(month|:,year|:)) * factor_rad
forrest@0
   195
forrest@0
   196
    data = fo->LE_f
forrest@0
   197
    data_ob_ann(2,:,n) = dim_avg(data(month|:,year|:))
forrest@0
   198
forrest@0
   199
    data = fo->H_f
forrest@0
   200
    data_ob_ann(3,:,n) = dim_avg(data(month|:,year|:))
forrest@0
   201
forrest@0
   202
    data = fo->GPP_or_MDS
forrest@0
   203
    data_ob_ann(4,:,n) = dim_avg(data(month|:,year|:))
forrest@0
   204
forrest@0
   205
    data = fo->Reco_or
forrest@0
   206
    data_ob_ann(5,:,n) = dim_avg(data(month|:,year|:))
forrest@0
   207
forrest@0
   208
    delete (data) 
forrest@0
   209
    delete (fo)   
forrest@0
   210
 end do
forrest@0
   211
forrest@0
   212
;--------------------------------------------------------------
forrest@0
   213
; find (# of year observed) >=4 and year_ob_i <= 2001
forrest@0
   214
forrest@0
   215
 i_long_ob = ind(year_station .ge. 4 .and. year_ob_i .le. 2001)
forrest@0
   216
 
forrest@0
   217
 station_long      = station(i_long_ob)
forrest@0
   218
 lat_ob_long       = lat_ob(i_long_ob)
forrest@0
   219
 lon_ob_long       = lat_ob(i_long_ob)
forrest@0
   220
 year_ob_long      = year_ob(i_long_ob)
forrest@0
   221
 year_ob_i_long    = year_ob_i(i_long_ob)
forrest@0
   222
 year_ob_f_long    = year_ob_f(i_long_ob)
forrest@0
   223
 year_station_long = year_station(i_long_ob)
forrest@0
   224
forrest@0
   225
 nstation_long     = dimsizes(station_long)
forrest@0
   226
forrest@0
   227
;=========================================================
forrest@0
   228
;   get model data at observed lat-lon
forrest@0
   229
forrest@0
   230
    fm   = addfile (dirm+film8,"r")
forrest@0
   231
forrest@0
   232
    xm   = fm->lon
forrest@0
   233
    ym   = fm->lat 
forrest@0
   234
    date = fm->date
forrest@0
   235
forrest@0
   236
    date_dim = dimsizes(date)
forrest@0
   237
    nyear  = date_dim(0)
forrest@0
   238
forrest@0
   239
    data_mod      = new ((/nfield,nyear,nmonth,nstation/),float)
forrest@0
   240
    data_mod_ann  = new ((/nfield,nmonth,nstation/),float)
forrest@0
   241
    data_mod_long = new ((/nfield,nyear,nmonth,nstation_long/),float)
forrest@0
   242
forrest@0
   243
;   change to unit of observed (u mol/m2/s)
forrest@0
   244
;   Model_units [=] gC/m2/s
forrest@0
   245
;   12. = molecular weight of C
forrest@0
   246
;   u mol = 1e-6 mol
forrest@0
   247
forrest@0
   248
    factor = 1.e6 /12.
forrest@0
   249
forrest@0
   250
;------------------------------------------------------------
forrest@0
   251
;   interpolate model data into observed station
forrest@0
   252
;   note: model is 0-360E, 90S-90N
forrest@0
   253
forrest@0
   254
;   to be able to handle observation at (-89.98,-24.80)
forrest@0
   255
    ym(0) = -90.
forrest@0
   256
;------------------------------------------------------------  
forrest@0
   257
forrest@0
   258
if (ENERGY .eq. "old") then
forrest@0
   259
forrest@0
   260
  data = fm->NEE
forrest@0
   261
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   262
  data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux
forrest@0
   263
  data_mod(0,:,:,:) = yy(:,:,:) * factor_flux
forrest@0
   264
forrest@0
   265
;;data  = fm->NETRAD
forrest@0
   266
; data  = fm->FSA
forrest@0
   267
; data1 = fm->FIRA
forrest@0
   268
; data  = data - data1
forrest@0
   269
; delete (data1)
forrest@0
   270
forrest@0
   271
  data  = fm->FSDS
forrest@0
   272
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   273
  data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:))
forrest@0
   274
  data_mod(1,:,:,:) = yy(:,:,:)  
forrest@0
   275
forrest@0
   276
forrest@0
   277
; data  = fm->LATENT
forrest@0
   278
  data  = fm->FCEV
forrest@0
   279
  data1 = fm->FCTR
forrest@0
   280
  data2 = fm->FGEV
forrest@0
   281
  data  = data + data1 + data2
forrest@0
   282
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   283
  data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:))
forrest@0
   284
  data_mod(2,:,:,:) = yy(:,:,:) 
forrest@0
   285
  delete (data1)
forrest@0
   286
  delete (data2)
forrest@0
   287
 
forrest@0
   288
; data = fm->SENSIBLE
forrest@0
   289
  data  = fm->FSH
forrest@0
   290
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   291
  data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 
forrest@0
   292
  data_mod(3,:,:,:) = yy(:,:,:)  
forrest@0
   293
forrest@0
   294
else
forrest@0
   295
forrest@0
   296
  data = fm->NEE
forrest@0
   297
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   298
  data_mod_ann(0,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux
forrest@0
   299
  data_mod(0,:,:,:) = yy(:,:,:) * factor_flux
forrest@0
   300
forrest@0
   301
; data = fm->NETRAD
forrest@0
   302
  data = fm->FSDS
forrest@0
   303
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   304
  data_mod_ann(1,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 
forrest@0
   305
  data_mod(1,:,:,:) = yy(:,:,:) 
forrest@0
   306
forrest@0
   307
  data = fm->LATENT
forrest@0
   308
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   309
  data_mod_ann(2,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 
forrest@0
   310
  data_mod(2,:,:,:) = yy(:,:,:) 
forrest@0
   311
forrest@0
   312
; data = fm->SENSIBLE
forrest@0
   313
  data = fm->FSH
forrest@0
   314
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   315
  data_mod_ann(3,:,:)= dim_avg(yy(month|:,pts|:,year|:)) 
forrest@0
   316
  data_mod(3,:,:,:) = yy(:,:,:) 
forrest@0
   317
forrest@0
   318
end if
forrest@0
   319
forrest@0
   320
  data = fm->GPP
forrest@0
   321
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   322
  data_mod_ann(4,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux
forrest@0
   323
  data_mod(4,:,:,:) = yy(:,:,:) * factor_flux 
forrest@0
   324
forrest@0
   325
  data = fm->ER
forrest@0
   326
  yy = linint2_points_Wrap(xm,ym,data,True,lon_ob,lat_ob,0)
forrest@0
   327
  data_mod_ann(5,:,:)= dim_avg(yy(month|:,pts|:,year|:)) * factor_flux
forrest@0
   328
  data_mod(5,:,:,:) = yy(:,:,:) * factor_flux
forrest@0
   329
forrest@0
   330
  data_mod_long(:,:,:,:) = data_mod(:,:,:,i_long_ob)
forrest@0
   331
forrest@0
   332
  delete (data_mod)
forrest@0
   333
  delete (fm) 
forrest@0
   334
  delete (data)
forrest@0
   335
  delete (yy)
forrest@0
   336
forrest@0
   337
  asciiwrite(station(0)+"_ob.txt", data_ob_ann(:,:,0))
forrest@0
   338
  asciiwrite(station(0)+"_"+model_name+".txt", data_mod_ann(:,:,0))
forrest@0
   339
forrest@0
   340
  asciiwrite(station(7)+"_ob.txt", data_ob_ann(:,:,7))
forrest@0
   341
  asciiwrite(station(7)+"_"+model_name+".txt", data_mod_ann(:,:,7))
forrest@0
   342
forrest@0
   343
  asciiwrite(station(21)+"_ob.txt", data_ob_ann(:,:,21))
forrest@0
   344
  asciiwrite(station(21)+"_"+model_name+".txt", data_mod_ann(:,:,21))
forrest@0
   345
forrest@0
   346
  asciiwrite(station(44)+"_ob.txt", data_ob_ann(:,:,44))
forrest@0
   347
  asciiwrite(station(44)+"_"+model_name+".txt", data_mod_ann(:,:,44))
forrest@0
   348
forrest@0
   349
  asciiwrite(station(50)+"_ob.txt", data_ob_ann(:,:,50))
forrest@0
   350
  asciiwrite(station(50)+"_"+model_name+".txt", data_mod_ann(:,:,50))
forrest@0
   351
forrest@0
   352
  asciiwrite(station(54)+"_ob.txt", data_ob_ann(:,:,54))
forrest@0
   353
  asciiwrite(station(54)+"_"+model_name+".txt", data_mod_ann(:,:,54))
forrest@0
   354
forrest@0
   355
  asciiwrite(station(62)+"_ob.txt", data_ob_ann(:,:,62))
forrest@0
   356
  asciiwrite(station(62)+"_"+model_name+".txt", data_mod_ann(:,:,62))
forrest@0
   357
exit
forrest@0
   358
forrest@0
   359
;*******************************************************************
forrest@0
   360
; for station line plot
forrest@0
   361
;*******************************************************************
forrest@0
   362
forrest@0
   363
; for x-axis in xyplot
forrest@0
   364
  mon = ispan(1,12,1)
forrest@0
   365
  mon@long_name = "month"
forrest@0
   366
forrest@0
   367
  res                   = True               ; plot mods desired
forrest@0
   368
  res@xyLineThicknesses = (/2.0,2.0/)        ; make 2nd lines thicker
forrest@0
   369
  res@xyLineColors      = (/"blue","red"/)   ; line color (ob,model)
forrest@0
   370
forrest@0
   371
  res@tmXBFormat  = "f"                      ; not to add trailing zeros
forrest@0
   372
forrest@0
   373
;-------------------------------------------------------------------------
forrest@0
   374
; Add a boxed legend using the more simple method
forrest@0
   375
forrest@0
   376
  res@pmLegendDisplayMode    = "Always"
forrest@0
   377
; res@pmLegendWidthF         = 0.1
forrest@0
   378
  res@pmLegendWidthF         = 0.08
forrest@0
   379
  res@pmLegendHeightF        = 0.06
forrest@0
   380
; res@pmLegendOrthogonalPosF = -1.17
forrest@0
   381
; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
forrest@0
   382
  res@pmLegendOrthogonalPosF = -0.30  ;(downward)
forrest@0
   383
forrest@0
   384
; res@pmLegendParallelPosF   =  0.18
forrest@0
   385
  res@pmLegendParallelPosF   =  0.23  ;(rightward)
forrest@0
   386
forrest@0
   387
; res@lgPerimOn             = False
forrest@0
   388
  res@lgLabelFontHeightF     = 0.015
forrest@0
   389
  res@xyExplicitLegendLabels = (/"observed",model_name/)
forrest@0
   390
;-------------------------------------------------------------------
forrest@0
   391
; for panel plot
forrest@0
   392
  res@gsnFrame     = False                   ; Do not draw plot 
forrest@0
   393
  res@gsnDraw      = False                   ; Do not advance frame
forrest@0
   394
forrest@0
   395
  pres                            = True     ; panel plot mods desired
forrest@0
   396
  pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
forrest@0
   397
                                             ; indiv. plots in panel
forrest@0
   398
  pres@gsnMaximize                = True     ; fill the page
forrest@0
   399
;-------------------------------------------------------------------
forrest@0
   400
forrest@0
   401
;==============================================================
forrest@0
   402
; get ob data at each site with long observation
forrest@0
   403
forrest@0
   404
  do n = 0,nstation_long-1
forrest@0
   405
forrest@0
   406
;##################################################################
forrest@0
   407
; hardwired: model    up to year 2004
forrest@0
   408
;            observed up to year 2006
forrest@0
   409
forrest@0
   410
    year_setback = 0
forrest@0
   411
forrest@0
   412
    nyear = year_station_long(n)
forrest@0
   413
forrest@0
   414
    if (year_ob_f_long(n).eq. 2006) then
forrest@0
   415
       year_setback = 2006 -2004
forrest@0
   416
    end if
forrest@0
   417
    if (year_ob_f_long(n).eq. 2005) then
forrest@0
   418
       year_setback = 2005 -2004
forrest@0
   419
    end if
forrest@0
   420
;##################################################################
forrest@0
   421
    
forrest@0
   422
    ntime = (nyear - year_setback) * nmonth
forrest@0
   423
forrest@0
   424
    data_ob   = new ((/nfield, nyear, nmonth/),float)
forrest@0
   425
forrest@0
   426
    dir_f = dir_root + station_long(n)+"/"
forrest@0
   427
    fil_f = "timeseries_L4_m.nc"
forrest@0
   428
    fo    = addfile (dir_f+fil_f,"r")
forrest@0
   429
forrest@0
   430
    data_ob(0,:,:) = fo->NEE_or_fMDS
forrest@0
   431
    data_ob(1,:,:) = fo->Rg_f
forrest@0
   432
    data_ob(2,:,:) = fo->LE_f
forrest@0
   433
    data_ob(3,:,:) = fo->H_f
forrest@0
   434
    data_ob(4,:,:) = fo->GPP_or_MDS
forrest@0
   435
    data_ob(5,:,:) = fo->Reco_or
forrest@0
   436
forrest@0
   437
    data_ob(1,:,:) = data_ob(1,:,:) * factor_rad
forrest@0
   438
forrest@0
   439
    delete (fo)
forrest@0
   440
forrest@0
   441
    timeI = new((/ntime/),integer)
forrest@0
   442
    timeF = new((/ntime/),float)
forrest@0
   443
    timeI = ispan(1,ntime,1)
forrest@0
   444
    timeF = year_ob_i_long(n) + (timeI-1)/12.
forrest@0
   445
    timeF@long_name = "year" 
forrest@0
   446
forrest@0
   447
    plot_data = new((/2,ntime/),float)
forrest@0
   448
 
forrest@0
   449
;----------------------------
forrest@0
   450
; for model_vs_ob
forrest@0
   451
forrest@0
   452
    plot_name = station_long(n)+"_tseries_vs_ob"
forrest@0
   453
    title = station_long(n)+"("+sprintf("%5.2f",lat_ob_long(n))+","+sprintf("%5.2f",lon_ob_long(n))+")"    
forrest@0
   454
    res@tiMainString = title
forrest@0
   455
forrest@0
   456
    wks = gsn_open_wks (plot_type,plot_name)
forrest@0
   457
    plot=new(nfield,graphic)                         ; create graphic array   
forrest@0
   458
forrest@0
   459
    i_year_mod_i = year_ob_i_long(n) - 1990 
forrest@0
   460
    i_year_mod_f = i_year_mod_i + nyear - 1 - year_setback
forrest@0
   461
forrest@0
   462
    i_year_ob_f =  nyear - year_setback - 1 
forrest@0
   463
forrest@0
   464
;   print (nyear)
forrest@0
   465
;   print (i_year_ob_f)
forrest@0
   466
;   print (i_year_mod_i)
forrest@0
   467
;   print (i_year_mod_f)
forrest@0
   468
forrest@0
   469
    do i = 0,nfield-1                           
forrest@0
   470
       plot_data(0,:) = ndtooned(data_ob (i,0:i_year_ob_f,:))
forrest@0
   471
       plot_data(1,:) = ndtooned(data_mod_long(i,i_year_mod_i:i_year_mod_f,:,n))
forrest@0
   472
       plot_data@long_name = field(i)+" ("+field_unit(i)+")"   
forrest@0
   473
       plot(i)=gsn_csm_xy(wks,timeF,plot_data,res)            ; create plot 
forrest@0
   474
    end do
forrest@0
   475
   
forrest@0
   476
    gsn_panel(wks,plot,(/3,2/),pres)                    ; create panel plot
forrest@0
   477
forrest@0
   478
    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   479
           "rm "+plot_name+"."+plot_type)
forrest@0
   480
forrest@0
   481
    clear (wks)  
forrest@0
   482
    delete (plot)
forrest@0
   483
forrest@0
   484
    delete (data_ob)
forrest@0
   485
    delete (timeI)    
forrest@0
   486
    delete (timeF)
forrest@0
   487
    delete (plot_data)
forrest@0
   488
forrest@0
   489
 end do
forrest@0
   490
forrest@0
   491
;###################################################################
forrest@0
   492
; for the following tables,
forrest@0
   493
; sort by latitude in decending order (N->S)
forrest@0
   494
forrest@0
   495
  isort = dim_pqsort(lat_ob_long,-1)
forrest@0
   496
forrest@0
   497
  station_sort = station_long(isort)
forrest@0
   498
  year_ob_sort = year_ob_long(isort)
forrest@0
   499
  lat_ob_sort  = lat_ob_long(isort)
forrest@0
   500
  lon_ob_sort  = lon_ob_long(isort)
forrest@0
   501
 
forrest@0
   502
; print(isort)
forrest@0
   503
; print(lat_ob_sort)
forrest@0
   504
forrest@0
   505
;*******************************************************************
forrest@0
   506
; html table of site: observed
forrest@0
   507
;*******************************************************************
forrest@0
   508
  output_html = "tseries_vs_ob.html"
forrest@0
   509
forrest@0
   510
  header = (/"<HTML>" \
forrest@0
   511
            ,"<HEAD>" \
forrest@0
   512
            ,"<TITLE>CLAMP metrics</TITLE>" \
forrest@0
   513
            ,"</HEAD>" \
forrest@0
   514
            ,"<H1>Timeseries at Site: "+model_name+" vs Observation</H1>" \
forrest@0
   515
            /) 
forrest@0
   516
  footer = "</HTML>"
forrest@0
   517
forrest@0
   518
  table_header = (/ \
forrest@0
   519
        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
forrest@0
   520
       ,"<tr>" \
forrest@0
   521
       ,"   <th bgcolor=DDDDDD >Site Name</th>" \
forrest@0
   522
       ,"   <th bgcolor=DDDDDD >Latitude</th>" \
forrest@0
   523
       ,"   <th bgcolor=DDDDDD >Longitude</th>" \
forrest@0
   524
       ,"   <th bgcolor=DDDDDD >Observed</th>" \ 
forrest@0
   525
       ,"</tr>" \
forrest@0
   526
       /)
forrest@0
   527
  table_footer = "</table>"
forrest@0
   528
  row_header = "<tr>"
forrest@0
   529
  row_footer = "</tr>"
forrest@0
   530
forrest@0
   531
  lines = new(50000,string)
forrest@0
   532
  nline = 0
forrest@0
   533
forrest@0
   534
  set_line(lines,nline,header)
forrest@0
   535
  set_line(lines,nline,table_header)
forrest@0
   536
;-----------------------------------------------
forrest@0
   537
; row of table
forrest@0
   538
  
forrest@0
   539
  do n = 0,nstation_long-1
forrest@0
   540
forrest@0
   541
     set_line(lines,nline,row_header)
forrest@0
   542
forrest@0
   543
     txt0 = station_sort(n)
forrest@0
   544
     txt1 = sprintf("%5.2f", lat_ob_sort(n))
forrest@0
   545
     txt2 = sprintf("%5.2f", lon_ob_sort(n))
forrest@0
   546
     txt3 = year_ob_sort(n)
forrest@0
   547
forrest@0
   548
     set_line(lines,nline,"<th><a href="+txt0+"_tseries_vs_ob.png>"+txt0+"</a></th>")
forrest@0
   549
     set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0
   550
     set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0
   551
     set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0
   552
forrest@0
   553
     set_line(lines,nline,row_footer)
forrest@0
   554
  end do
forrest@0
   555
;-----------------------------------------------
forrest@0
   556
  set_line(lines,nline,table_footer)
forrest@0
   557
  set_line(lines,nline,footer) 
forrest@0
   558
forrest@0
   559
; Now write to an HTML file.
forrest@0
   560
  idx = ind(.not.ismissing(lines))
forrest@0
   561
  if(.not.any(ismissing(idx))) then
forrest@0
   562
    asciiwrite(output_html,lines(idx))
forrest@0
   563
  else
forrest@0
   564
   print ("error?")
forrest@0
   565
  end if
forrest@0
   566
  delete (idx)
forrest@0
   567
forrest@0
   568
  delete (isort)
forrest@0
   569
  delete (station_sort)
forrest@0
   570
  delete (year_ob_sort)
forrest@0
   571
  delete (lat_ob_sort)
forrest@0
   572
  delete (lon_ob_sort)
forrest@0
   573
forrest@0
   574
;************************************************************
forrest@0
   575
; compute annual cycle correlation coef and M score
forrest@0
   576
;************************************************************
forrest@0
   577
forrest@0
   578
 score_max = 1.
forrest@0
   579
forrest@0
   580
 ccr     = new ((/nstation, nfield/),float)
forrest@0
   581
 M_score = new ((/nstation, nfield/),float) 
forrest@0
   582
forrest@0
   583
 do n=0,nstation-1
forrest@0
   584
 do m=0,nfield-1   
forrest@0
   585
    ccr(n,m) = esccr(data_ob_ann(m,:,n),data_mod_ann(m,:,n),0)
forrest@0
   586
    bias = sum(abs(data_mod_ann(m,:,n)-data_ob_ann(m,:,n))/(abs(data_mod_ann(m,:,n))+abs(data_ob_ann(m,:,n))))
forrest@0
   587
    M_score(n,m) = (1. -(bias/nmonth)) * score_max
forrest@0
   588
 end do
forrest@0
   589
 end do
forrest@0
   590
forrest@0
   591
 M_nee = avg(M_score(:,0))
forrest@0
   592
 M_rad = avg(M_score(:,1))
forrest@0
   593
 M_lh  = avg(M_score(:,2))
forrest@0
   594
 M_sh  = avg(M_score(:,3))
forrest@0
   595
 M_gpp = avg(M_score(:,4))
forrest@0
   596
 M_er  = avg(M_score(:,5))
forrest@0
   597
 M_all = M_nee+ M_rad +M_lh + M_sh + M_gpp + M_er
forrest@0
   598
forrest@0
   599
 M_ameriflux_nee = sprintf("%.2f", M_nee)
forrest@0
   600
 M_ameriflux_rad = sprintf("%.2f", M_rad)
forrest@0
   601
 M_ameriflux_lh  = sprintf("%.2f", M_lh )
forrest@0
   602
 M_ameriflux_sh  = sprintf("%.2f", M_sh )
forrest@0
   603
 M_ameriflux_gpp = sprintf("%.2f", M_gpp)
forrest@0
   604
 M_ameriflux_er  = sprintf("%.2f", M_er )
forrest@0
   605
 M_ameriflux_all = sprintf("%.2f", M_all)
forrest@0
   606
forrest@0
   607
;*******************************************************************
forrest@0
   608
; for station line plot
forrest@0
   609
;*******************************************************************
forrest@0
   610
forrest@0
   611
; for x-axis in xyplot
forrest@0
   612
  mon = ispan(1,12,1)
forrest@0
   613
  mon@long_name = "month"
forrest@0
   614
forrest@0
   615
;-------------------------------------------------------------------
forrest@0
   616
forrest@0
   617
  plot_data   = new((/2,nmonth/),float)
forrest@0
   618
  plot_data!0 = "case"
forrest@0
   619
  plot_data!1 = "month"
forrest@0
   620
forrest@0
   621
  do n = 0,nstation-1
forrest@0
   622
;----------------------------
forrest@0
   623
; for observed
forrest@0
   624
forrest@0
   625
    plot_name = station(n)+"_ob"    
forrest@0
   626
    title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"
forrest@0
   627
    res@tiMainString = title
forrest@0
   628
forrest@0
   629
    wks = gsn_open_wks (plot_type,plot_name)
forrest@0
   630
    plot=new(nfield,graphic)                        ; create graphic array   
forrest@0
   631
forrest@0
   632
    do i = 0,nfield-1                           
forrest@0
   633
       plot_data(0,:) = (/data_ob_ann(i,:,n)/)
forrest@0
   634
       plot_data@long_name = field(i)+" ("+field_unit(i)+")"   
forrest@0
   635
       plot(i)=gsn_csm_xy(wks,mon,plot_data(0,:),res)    ; create plot 
forrest@0
   636
    end do
forrest@0
   637
forrest@0
   638
    gsn_panel(wks,plot,(/3,2/),pres)                 ; create panel plot
forrest@0
   639
forrest@0
   640
    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   641
           "rm "+plot_name+"."+plot_type)
forrest@0
   642
forrest@0
   643
    clear (wks)  
forrest@0
   644
    delete (plot)
forrest@0
   645
;----------------------------
forrest@0
   646
; for model_vs_ob
forrest@0
   647
forrest@0
   648
    plot_name = station(n)+"_model_vs_ob"
forrest@0
   649
    title = station(n)+"("+sprintf("%5.2f",lat_ob(n))+","+sprintf("%5.2f",lon_ob(n))+")"    
forrest@0
   650
    res@tiMainString = title
forrest@0
   651
forrest@0
   652
    wks = gsn_open_wks (plot_type,plot_name)
forrest@0
   653
    plot=new(nfield,graphic)                         ; create graphic array   
forrest@0
   654
forrest@0
   655
    do i = 0,nfield-1                           
forrest@0
   656
       plot_data(0,:) = (/data_ob_ann(i,:,n)/)
forrest@0
   657
       plot_data(1,:) = (/data_mod_ann(i,:,n)/)
forrest@0
   658
       plot_data@long_name = field(i)+" ("+field_unit(i)+")"   
forrest@0
   659
       plot(i)=gsn_csm_xy(wks,mon,plot_data,res)     ; create plot 
forrest@0
   660
    end do
forrest@0
   661
   
forrest@0
   662
    gsn_panel(wks,plot,(/3,2/),pres)                 ; create panel plot
forrest@0
   663
forrest@0
   664
    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   665
           "rm "+plot_name+"."+plot_type)
forrest@0
   666
forrest@0
   667
    clear (wks)  
forrest@0
   668
    delete (plot)
forrest@0
   669
 end do
forrest@0
   670
forrest@0
   671
;###################################################################
forrest@0
   672
; for the following tables,
forrest@0
   673
; sort by latitude in decending order (N->S)
forrest@0
   674
; sort by lat in decending order (N->S)
forrest@0
   675
forrest@0
   676
  isort = dim_pqsort(lat_ob,-1)
forrest@0
   677
forrest@0
   678
  station_sort = station(isort)
forrest@0
   679
  year_ob_sort = year_ob(isort)
forrest@0
   680
  lat_ob_sort  = lat_ob(isort)
forrest@0
   681
  lon_ob_sort  = lon_ob(isort)
forrest@0
   682
  M_score_sort = M_score(isort,:)
forrest@0
   683
 
forrest@0
   684
; print(isort)
forrest@0
   685
; print(lat_ob_sort)
forrest@0
   686
;###################################################################
forrest@0
   687
;*******************************************************************
forrest@0
   688
; html table of site: observed
forrest@0
   689
;*******************************************************************
forrest@0
   690
  output_html = "line_ob.html"
forrest@0
   691
forrest@0
   692
  header = (/"<HTML>" \
forrest@0
   693
            ,"<HEAD>" \
forrest@0
   694
            ,"<TITLE>CLAMP metrics</TITLE>" \
forrest@0
   695
            ,"</HEAD>" \
forrest@0
   696
            ,"<H1>Energy at Site: Observation</H1>" \
forrest@0
   697
            /) 
forrest@0
   698
  footer = "</HTML>"
forrest@0
   699
forrest@0
   700
  table_header = (/ \
forrest@0
   701
        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
forrest@0
   702
       ,"<tr>" \
forrest@0
   703
       ,"   <th bgcolor=DDDDDD >Site Name</th>" \
forrest@0
   704
       ,"   <th bgcolor=DDDDDD >Latitude</th>" \
forrest@0
   705
       ,"   <th bgcolor=DDDDDD >Longitude</th>" \
forrest@0
   706
       ,"   <th bgcolor=DDDDDD >Observed</th>" \ 
forrest@0
   707
       ,"</tr>" \
forrest@0
   708
       /)
forrest@0
   709
  table_footer = "</table>"
forrest@0
   710
  row_header = "<tr>"
forrest@0
   711
  row_footer = "</tr>"
forrest@0
   712
forrest@0
   713
  lines = new(50000,string)
forrest@0
   714
  nline = 0
forrest@0
   715
forrest@0
   716
  set_line(lines,nline,header)
forrest@0
   717
  set_line(lines,nline,table_header)
forrest@0
   718
;-----------------------------------------------
forrest@0
   719
; row of table
forrest@0
   720
  
forrest@0
   721
  do n = 0,nstation-1
forrest@0
   722
     set_line(lines,nline,row_header)
forrest@0
   723
forrest@0
   724
     txt0 = station_sort(n)
forrest@0
   725
     txt1 = sprintf("%5.2f", lat_ob_sort(n))
forrest@0
   726
     txt2 = sprintf("%5.2f", lon_ob_sort(n))
forrest@0
   727
     txt3 = year_ob_sort(n)
forrest@0
   728
forrest@0
   729
     set_line(lines,nline,"<th><a href="+txt0+"_ob.png>"+txt0+"</a></th>")
forrest@0
   730
     set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0
   731
     set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0
   732
     set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0
   733
forrest@0
   734
     set_line(lines,nline,row_footer)
forrest@0
   735
  end do
forrest@0
   736
;-----------------------------------------------
forrest@0
   737
  set_line(lines,nline,table_footer)
forrest@0
   738
  set_line(lines,nline,footer) 
forrest@0
   739
forrest@0
   740
; Now write to an HTML file.
forrest@0
   741
  idx = ind(.not.ismissing(lines))
forrest@0
   742
  if(.not.any(ismissing(idx))) then
forrest@0
   743
    asciiwrite(output_html,lines(idx))
forrest@0
   744
  else
forrest@0
   745
   print ("error?")
forrest@0
   746
  end if
forrest@0
   747
  delete (idx)
forrest@0
   748
forrest@0
   749
;*******************************************************************
forrest@0
   750
; score and line table : model vs observed
forrest@0
   751
;*******************************************************************
forrest@0
   752
  output_html = "score+line_vs_ob.html"
forrest@0
   753
forrest@0
   754
  header = (/"<HTML>" \
forrest@0
   755
            ,"<HEAD>" \
forrest@0
   756
            ,"<TITLE>CLAMP metrics</TITLE>" \
forrest@0
   757
            ,"</HEAD>" \
forrest@0
   758
            ,"<H1>Energy at Site: Model "+model_name+"</H1>" \
forrest@0
   759
            /) 
forrest@0
   760
  footer = "</HTML>"
forrest@0
   761
forrest@0
   762
  delete (table_header)
forrest@0
   763
  table_header = (/ \
forrest@0
   764
        "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
forrest@0
   765
       ,"<tr>" \
forrest@0
   766
       ,"   <th bgcolor=DDDDDD >Site Name</th>" \
forrest@0
   767
       ,"   <th bgcolor=DDDDDD >Latitude</th>" \
forrest@0
   768
       ,"   <th bgcolor=DDDDDD >Longitude</th>" \
forrest@0
   769
       ,"   <th bgcolor=DDDDDD >Observed</th>" \
forrest@0
   770
       ,"   <th bgcolor=DDDDDD >NEE Flux</th>" \
forrest@0
   771
       ,"   <th bgcolor=DDDDDD >Shortwave <br> Incoming</th>" \
forrest@0
   772
       ,"   <th bgcolor=DDDDDD >Latent Heat</th>" \
forrest@0
   773
       ,"   <th bgcolor=DDDDDD >Sensible Heat</th>" \
forrest@0
   774
       ,"   <th bgcolor=DDDDDD >GPP Flux</th>" \
forrest@0
   775
       ,"   <th bgcolor=DDDDDD >Respiration</th>" \
forrest@0
   776
       ,"   <th bgcolor=DDDDDD >Average</th>" \
forrest@0
   777
       ,"</tr>" \
forrest@0
   778
       /)
forrest@0
   779
  table_footer = "</table>"
forrest@0
   780
  row_header = "<tr>"
forrest@0
   781
  row_footer = "</tr>"
forrest@0
   782
forrest@0
   783
  lines = new(50000,string)
forrest@0
   784
  nline = 0
forrest@0
   785
forrest@0
   786
  set_line(lines,nline,header)
forrest@0
   787
  set_line(lines,nline,table_header)
forrest@0
   788
;-----------------------------------------------
forrest@0
   789
; row of table
forrest@0
   790
  
forrest@0
   791
  do n = 0,nstation-1
forrest@0
   792
     set_line(lines,nline,row_header)
forrest@0
   793
forrest@0
   794
     txt0  = station_sort(n)
forrest@0
   795
     txt1  = sprintf("%5.2f", lat_ob_sort(n))
forrest@0
   796
     txt2  = sprintf("%5.2f", lon_ob_sort(n))
forrest@0
   797
     txt3  = year_ob_sort(n)
forrest@0
   798
     txt4  = sprintf("%5.2f", M_score_sort(n,0))
forrest@0
   799
     txt5  = sprintf("%5.2f", M_score_sort(n,1))
forrest@0
   800
     txt6  = sprintf("%5.2f", M_score_sort(n,2))
forrest@0
   801
     txt7  = sprintf("%5.2f", M_score_sort(n,3))
forrest@0
   802
     txt8  = sprintf("%5.2f", M_score_sort(n,4))
forrest@0
   803
     txt9  = sprintf("%5.2f", M_score_sort(n,5))
forrest@0
   804
     txt10 = sprintf("%5.2f", avg(M_score_sort(n,:)))
forrest@0
   805
forrest@0
   806
     set_line(lines,nline,"<th><a href="+txt0+"_model_vs_ob.png>"+txt0+"</a></th>")
forrest@0
   807
     set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0
   808
     set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0
   809
     set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0
   810
     set_line(lines,nline,"<th>"+txt4+"</th>")
forrest@0
   811
     set_line(lines,nline,"<th>"+txt5+"</th>")
forrest@0
   812
     set_line(lines,nline,"<th>"+txt6+"</th>")
forrest@0
   813
     set_line(lines,nline,"<th>"+txt7+"</th>")
forrest@0
   814
     set_line(lines,nline,"<th>"+txt8+"</th>")
forrest@0
   815
     set_line(lines,nline,"<th>"+txt9+"</th>")
forrest@0
   816
     set_line(lines,nline,"<th>"+txt10+"</th>")
forrest@0
   817
forrest@0
   818
     set_line(lines,nline,row_footer)
forrest@0
   819
  end do
forrest@0
   820
forrest@0
   821
; last row, summary
forrest@0
   822
  set_line(lines,nline,row_header)
forrest@0
   823
forrest@0
   824
  txt0  = "All_"+sprintf("%.0f", nstation)
forrest@0
   825
  txt1  = "-"
forrest@0
   826
  txt2  = "-"
forrest@0
   827
  txt3  = "-"
forrest@0
   828
  txt4  = M_ameriflux_nee
forrest@0
   829
  txt5  = M_ameriflux_rad
forrest@0
   830
  txt6  = M_ameriflux_lh
forrest@0
   831
  txt7  = M_ameriflux_sh
forrest@0
   832
  txt8  = M_ameriflux_gpp
forrest@0
   833
  txt9  = M_ameriflux_er
forrest@0
   834
  txt10 = M_ameriflux_all
forrest@0
   835
forrest@0
   836
  set_line(lines,nline,"<th>"+txt0+"</th>")
forrest@0
   837
  set_line(lines,nline,"<th>"+txt1+"</th>")
forrest@0
   838
  set_line(lines,nline,"<th>"+txt2+"</th>")
forrest@0
   839
  set_line(lines,nline,"<th>"+txt3+"</th>")
forrest@0
   840
  set_line(lines,nline,"<th>"+txt4+"</th>")
forrest@0
   841
  set_line(lines,nline,"<th>"+txt5+"</th>")
forrest@0
   842
  set_line(lines,nline,"<th>"+txt6+"</th>")
forrest@0
   843
  set_line(lines,nline,"<th>"+txt7+"</th>")
forrest@0
   844
  set_line(lines,nline,"<th>"+txt8+"</th>")
forrest@0
   845
  set_line(lines,nline,"<th>"+txt9+"</th>")
forrest@0
   846
  set_line(lines,nline,"<th>"+txt10+"</th>")
forrest@0
   847
forrest@0
   848
  set_line(lines,nline,row_footer)
forrest@0
   849
;-----------------------------------------------
forrest@0
   850
  set_line(lines,nline,table_footer)
forrest@0
   851
  set_line(lines,nline,footer) 
forrest@0
   852
forrest@0
   853
; Now write to an HTML file.
forrest@0
   854
  idx = ind(.not.ismissing(lines))
forrest@0
   855
  if(.not.any(ismissing(idx))) then
forrest@0
   856
    asciiwrite(output_html,lines(idx))
forrest@0
   857
  else
forrest@0
   858
   print ("error?")
forrest@0
   859
  end if
forrest@0
   860
  delete (idx)
forrest@0
   861
forrest@0
   862
;**************************************************************************************
forrest@0
   863
; update score
forrest@0
   864
;**************************************************************************************
forrest@0
   865
 
forrest@0
   866
  if (isvar("compare")) then
forrest@0
   867
     system("sed 1,/M_ameriflux_nee/s//"+M_ameriflux_nee+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0
   868
            "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0
   869
            "sed 1,/M_ameriflux_rad/s//"+M_ameriflux_rad+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0
   870
            "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0
   871
            "sed 1,/M_ameriflux_lh/s//"+M_ameriflux_lh+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0
   872
            "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0
   873
            "sed 1,/M_ameriflux_sh/s//"+M_ameriflux_sh+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0
   874
            "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0
   875
            "sed 1,/M_ameriflux_gpp/s//"+M_ameriflux_gpp+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0
   876
            "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0
   877
            "sed 1,/M_ameriflux_er/s//"+M_ameriflux_er+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0
   878
            "mv -f "+html_new2+" "+html_name2)
forrest@0
   879
  end if
forrest@0
   880
forrest@0
   881
  system("sed s#M_ameriflux_nee#"+M_ameriflux_nee+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   882
         "mv -f "+html_new+" "+html_name+";"+ \
forrest@0
   883
         "sed s#M_ameriflux_rad#"+M_ameriflux_rad+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   884
         "mv -f "+html_new+" "+html_name+";"+ \
forrest@0
   885
         "sed s#M_ameriflux_lh#"+M_ameriflux_lh+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   886
         "mv -f "+html_new+" "+html_name+";"+ \
forrest@0
   887
         "sed s#M_ameriflux_sh#"+M_ameriflux_sh+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   888
         "mv -f "+html_new+" "+html_name+";"+ \
forrest@0
   889
         "sed s#M_ameriflux_gpp#"+M_ameriflux_gpp+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   890
         "mv -f "+html_new+" "+html_name+";"+ \
forrest@0
   891
         "sed s#M_ameriflux_er#"+M_ameriflux_er+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   892
         "mv -f "+html_new+" "+html_name) 
forrest@0
   893
forrest@0
   894
;***************************************************************************
forrest@0
   895
; add total score and write to file
forrest@0
   896
;***************************************************************************
forrest@0
   897
  M_total = M_ameriflux_all
forrest@0
   898
forrest@0
   899
  asciiwrite("M_save.ameriflux", M_total)
forrest@0
   900
forrest@0
   901
;***************************************************************************
forrest@0
   902
; output plot and html
forrest@0
   903
;***************************************************************************
forrest@0
   904
  output_dir = model_name+"/ameriflux"
forrest@0
   905
forrest@0
   906
  system("mv *.png *.html " + output_dir) 
forrest@0
   907
;***************************************************************************
forrest@0
   908
forrest@0
   909
end
forrest@0
   910