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