all/04.biomass.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Thu, 26 Mar 2009 14:02:21 -0400
changeset 1 4be95183fbcd
parent 0 0c6405ab2ff4
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
begin
forrest@0
     7
forrest@0
     8
 plot_type     = "ps"
forrest@0
     9
 plot_type_new = "png"
forrest@0
    10
forrest@0
    11
;------------------------------------------------------
forrest@0
    12
; edit table.html of current model for movel1_vs_model2
forrest@0
    13
forrest@0
    14
 if (isvar("compare")) then
forrest@0
    15
    html_name2 = compare+"/table.html"  
forrest@0
    16
    html_new2  = html_name2 +".new"
forrest@0
    17
 end if
forrest@0
    18
forrest@0
    19
;------------------------------------------------------
forrest@0
    20
; edit table.html for current model
forrest@0
    21
forrest@0
    22
 html_name = model_name+"/table.html"  
forrest@0
    23
 html_new  = html_name +".new"
forrest@0
    24
forrest@0
    25
;---------------------------------
forrest@0
    26
; get model data: landfrac and area
forrest@0
    27
forrest@0
    28
  film_l   = "lnd_"+ model_grid + ".nc"
forrest@0
    29
  fm_l     = addfile (dirs+film_l,"r")  
forrest@0
    30
  
forrest@0
    31
  landfrac = fm_l->landfrac
forrest@0
    32
  area0    = fm_l->area
forrest@0
    33
forrest@0
    34
  delete (fm_l)
forrest@0
    35
forrest@0
    36
; change area from km**2 to m**2
forrest@0
    37
  area0 = area0 * 1.e6
forrest@0
    38
             
forrest@0
    39
;-----------------------------
forrest@0
    40
; take into account landfrac
forrest@0
    41
forrest@0
    42
  area0     = area0 * landfrac
forrest@0
    43
forrest@0
    44
;--------------------------------------------
forrest@0
    45
; read model data
forrest@0
    46
forrest@1
    47
 fm   = addfile (dirm+film11,"r")
forrest@0
    48
forrest@0
    49
 if (BGC .eq. "cn") then 
forrest@0
    50
    data1  = fm->LIVESTEMC
forrest@0
    51
    data2  = fm->DEADSTEMC
forrest@0
    52
    data3  = fm->LEAFC
forrest@0
    53
    datamod0 = data1(0,:,:)
forrest@0
    54
    datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
forrest@0
    55
 end if
forrest@0
    56
forrest@0
    57
 if (BGC .eq. "casa") then
forrest@0
    58
    factor_WOODC = 0.7  
forrest@0
    59
    data1  = fm->WOODC
forrest@0
    60
    data2  = fm->LEAFC
forrest@0
    61
    datamod0 = data1(0,:,:)
forrest@0
    62
    datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
forrest@0
    63
 end if
forrest@0
    64
forrest@0
    65
; unit: gC/m2
forrest@0
    66
forrest@0
    67
 xm       = fm->lon  
forrest@0
    68
 ym       = fm->lat
forrest@0
    69
forrest@0
    70
 delete (fm)
forrest@0
    71
forrest@0
    72
;------------------------------------------------
forrest@0
    73
; read amazon mask data
forrest@0
    74
forrest@0
    75
  dir_m = diro + "biomass/"
forrest@0
    76
  fil_m = "amazon_mask_"+ model_grid + ".nc"
forrest@0
    77
  fm    = addfile (dir_m+fil_m,"r")
forrest@0
    78
forrest@0
    79
  mask_amazon0 = fm->mask_amazon
forrest@0
    80
forrest@0
    81
  delete (fm)
forrest@0
    82
forrest@0
    83
;------------------------------------------------
forrest@0
    84
; read ob data
forrest@0
    85
forrest@1
    86
 ob_name = "LC15 Amazon Biomass"
forrest@0
    87
forrest@0
    88
 dir_b = diro + "biomass/"
forrest@0
    89
 fil_b = "amazon_biomass_"+model_grid+".nc"
forrest@0
    90
 fo   = addfile (dir_b+fil_b,"r")
forrest@0
    91
 
forrest@0
    92
 dataob   = fo->BIOMASS
forrest@0
    93
 xo       = fo->lon  
forrest@0
    94
 yo       = fo->lat
forrest@0
    95
forrest@0
    96
 delete (fo)
forrest@0
    97
forrest@0
    98
;************************************************
forrest@0
    99
; Units for these variables are:
forrest@0
   100
; dataob   : MgC/ha
forrest@0
   101
; datamod0 : gC/m2
forrest@1
   102
; We want to convert these to kgC/m2
forrest@0
   103
; ha = 100m*100m = 10,000 m2
forrest@1
   104
; MgC/ha*1000/10,000 = kgC/m2
forrest@0
   105
forrest@0
   106
  factor_aboveground = 0.5
forrest@0
   107
  factor_unit_ob     = 0.1
forrest@0
   108
  factor_unit_mod    = 0.001         
forrest@0
   109
forrest@0
   110
  dataob   = dataob * factor_aboveground * factor_unit_ob
forrest@0
   111
  datamod0 = datamod0 * factor_unit_mod 
forrest@0
   112
forrest@1
   113
  dataob@units      = "kg C m~S~-2~N~"
forrest@1
   114
  datamod0@units    = "kg C m~S~-2~N~"
forrest@0
   115
forrest@0
   116
  dataob@long_name      = "Amazon Biomass"
forrest@0
   117
  datamod0@long_name    = "Amazon Biomass"
forrest@0
   118
;********************************************************
forrest@0
   119
; get subset of datamod0 that match dataob
forrest@0
   120
  
forrest@0
   121
  nlon = dimsizes(xo)
forrest@0
   122
  nlat = dimsizes(yo)
forrest@0
   123
forrest@0
   124
  ind_lonL = ind(xm .eq. xo(0))
forrest@0
   125
  ind_lonR = ind(xm .eq. xo(nlon-1))
forrest@0
   126
  ind_latS = ind(ym .eq. yo(0))
forrest@0
   127
  ind_latN = ind(ym .eq. yo(nlat-1))
forrest@0
   128
forrest@0
   129
  datamod  = dataob
forrest@0
   130
  datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
forrest@0
   131
forrest@0
   132
  area  = dataob
forrest@0
   133
  area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR)
forrest@0
   134
forrest@0
   135
  mask_amazon  = dataob
forrest@0
   136
  mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR)
forrest@0
   137
forrest@0
   138
  mask_amazon@units = ""
forrest@0
   139
;********************************************************
forrest@0
   140
; sum over amazom_mask area:
forrest@0
   141
forrest@1
   142
; Peta g = 1.e15 g = 1.e12 kg
forrest@0
   143
  factor_unit = 1.e-12
forrest@0
   144
forrest@0
   145
; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
forrest@0
   146
forrest@0
   147
  Sum_area = sum(area*mask_amazon)*factor_unit
forrest@0
   148
forrest@0
   149
  Sum_ob  = sum(dataob*area*mask_amazon)*factor_unit
forrest@0
   150
  Sum_mod = sum(datamod*area*mask_amazon)*factor_unit
forrest@0
   151
forrest@0
   152
  avg_ob  = Sum_ob /Sum_area
forrest@0
   153
  avg_mod = Sum_mod/Sum_area
forrest@0
   154
forrest@0
   155
  Sum_biomass_ob  = sprintf("%.2f",Sum_ob )
forrest@0
   156
  Sum_biomass_mod = sprintf("%.2f",Sum_mod) 
forrest@0
   157
    
forrest@0
   158
;---------------------------------------------------------------------- 
forrest@0
   159
; contour plot res
forrest@0
   160
forrest@0
   161
  resg                     = True             ; Use plot options
forrest@0
   162
  resg@cnFillOn            = True             ; Turn on color fill
forrest@0
   163
  resg@gsnSpreadColors     = True             ; use full colormap
forrest@0
   164
  resg@cnLinesOn           = False            ; Turn off contourn lines
forrest@0
   165
  resg@mpFillOn            = False            ; Turn off map fill
forrest@0
   166
  resg@gsnAddCyclic        = False
forrest@0
   167
forrest@0
   168
  resg@gsnSpreadColors      = True            ; use full colormap
forrest@0
   169
  resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
forrest@0
   170
forrest@0
   171
  resg@mpMinLatF            = -21.1      ; range to zoom in on
forrest@0
   172
  resg@mpMaxLatF            =  13.8
forrest@0
   173
  resg@mpMinLonF            =  277.28
forrest@0
   174
  resg@mpMaxLonF            =  326.43
forrest@0
   175
;------------------------------------------------------------------------
forrest@0
   176
; mask plot
forrest@0
   177
forrest@0
   178
  plot_name = "mask_ob"
forrest@0
   179
forrest@0
   180
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
   181
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
   182
forrest@0
   183
;-----------------------------------------
forrest@0
   184
; plot area sum
forrest@0
   185
forrest@0
   186
  gRes  = True
forrest@0
   187
  gRes@txFontHeightF = 0.02
forrest@0
   188
; gRes@txAngleF = 90
forrest@0
   189
forrest@0
   190
  area_sum_text = "(mask area = "+sprintf("%.2f", Sum_area)+"e12 m2)"
forrest@0
   191
forrest@0
   192
  gsn_text_ndc(wks,area_sum_text,0.50,0.80,gRes)
forrest@0
   193
;-----------------------------------------
forrest@0
   194
forrest@0
   195
  resg@cnMinLevelValF      = 0.              ; Min level
forrest@0
   196
  resg@cnMaxLevelValF      = 1.              ; Max level
forrest@0
   197
  resg@cnLevelSpacingF     = 0.1             ; interval
forrest@0
   198
forrest@0
   199
  resg@tiMainString        = "Amazon Mask: grid = "+ model_grid
forrest@0
   200
  
forrest@0
   201
  plot = gsn_csm_contour_map_ce(wks,mask_amazon,resg)   
forrest@0
   202
forrest@0
   203
  delete (wks)
forrest@0
   204
forrest@0
   205
  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   206
         "rm "+plot_name+"."+plot_type)
forrest@1
   207
  ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0
   208
forrest@0
   209
;------------------------------------------------------------------------
forrest@0
   210
; contour ob
forrest@0
   211
forrest@0
   212
  resg@cnMinLevelValF       = 0.              ; Min level
forrest@0
   213
  resg@cnMaxLevelValF       = 30.             ; Max level
forrest@0
   214
  resg@cnLevelSpacingF      = 2.              ; interval
forrest@0
   215
  
forrest@0
   216
  plot_name = "global_ob"
forrest@0
   217
  title     = ob_name+" "+ model_grid
forrest@0
   218
  resg@tiMainString  = title
forrest@0
   219
forrest@0
   220
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
   221
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
   222
forrest@0
   223
  plot = gsn_csm_contour_map_ce(wks,dataob,resg)   
forrest@0
   224
forrest@0
   225
  delete (wks)
forrest@0
   226
forrest@0
   227
  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   228
         "rm "+plot_name+"."+plot_type)
forrest@1
   229
  ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0
   230
forrest@0
   231
;------------------------------------------------------------------------
forrest@0
   232
; contour model
forrest@0
   233
forrest@0
   234
  resg@cnMinLevelValF       = 0.              ; Min level
forrest@0
   235
  resg@cnMaxLevelValF       = 30.             ; Max level
forrest@0
   236
  resg@cnLevelSpacingF      = 2.              ; interval
forrest@0
   237
forrest@0
   238
  plot_name = "global_model"
forrest@1
   239
  title     = "Model "+ model_name + " (2000)"
forrest@0
   240
  resg@tiMainString  = title
forrest@0
   241
forrest@0
   242
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
   243
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
   244
forrest@0
   245
  plot = gsn_csm_contour_map_ce(wks,datamod,resg)   
forrest@0
   246
forrest@0
   247
  delete (wks)
forrest@0
   248
forrest@0
   249
  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   250
         "rm "+plot_name+"."+plot_type)
forrest@1
   251
  ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0
   252
forrest@0
   253
;------------------------------------------------------------------------
forrest@0
   254
; contour model vs ob
forrest@0
   255
forrest@0
   256
  plot_name = "global_model_vs_ob"
forrest@0
   257
forrest@0
   258
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
   259
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
   260
forrest@0
   261
  delete (plot)
forrest@0
   262
  plot=new(3,graphic)                        ; create graphic array
forrest@0
   263
forrest@0
   264
  resg@gsnFrame             = False          ; Do not draw plot 
forrest@0
   265
  resg@gsnDraw              = False          ; Do not advance frame
forrest@0
   266
forrest@0
   267
;(d) compute correlation coef and M score
forrest@0
   268
forrest@0
   269
  uu = ndtooned(datamod)
forrest@0
   270
  vv = ndtooned(dataob)
forrest@0
   271
 
forrest@0
   272
  good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))
forrest@0
   273
forrest@0
   274
  ug = uu(good)
forrest@0
   275
  vg = vv(good)
forrest@0
   276
forrest@0
   277
  ccrG = esccr(ug,vg,0)
forrest@0
   278
forrest@1
   279
  score_max = 10.0
forrest@0
   280
forrest@0
   281
; Miomass = (ccrG*ccrG)* score_max
forrest@0
   282
; new eq
forrest@0
   283
  bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
forrest@1
   284
  print ("bias=" + bias + "dimsizes(ug)=" + dimsizes(ug))
forrest@1
   285
  print ("other bias=" + sum(abs(ug-vg)/(ug+vg)))
forrest@0
   286
  Mbiomass  = (1. - (bias/dimsizes(ug)))*score_max
forrest@1
   287
  print ("Mbiomass=" + Mbiomass)
forrest@0
   288
  M_biomass = sprintf("%.2f", Mbiomass)
forrest@1
   289
  print ("M_biomass=" + M_biomass)
forrest@0
   290
forrest@0
   291
  if (isvar("compare")) then
forrest@0
   292
     system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \ 
forrest@0
   293
            "mv -f "+html_new2+" "+html_name2)
forrest@0
   294
  end if
forrest@0
   295
forrest@0
   296
  system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   297
         "mv -f "+html_new+" "+html_name)
forrest@0
   298
forrest@0
   299
; plot correlation coef
forrest@0
   300
forrest@0
   301
  gRes  = True
forrest@0
   302
  gRes@txFontHeightF = 0.02
forrest@0
   303
  gRes@txAngleF = 90
forrest@0
   304
forrest@0
   305
  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
forrest@0
   306
forrest@0
   307
  gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0
   308
;--------------------------------------------------------------------
forrest@0
   309
  
forrest@0
   310
;(a) ob
forrest@0
   311
forrest@0
   312
  title     = ob_name+" "+ model_grid
forrest@0
   313
  resg@tiMainString  = title
forrest@0
   314
forrest@0
   315
  plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)       
forrest@0
   316
forrest@0
   317
;(b) model
forrest@0
   318
forrest@1
   319
  title     = "Model "+ model_name + " (2000)"
forrest@0
   320
  resg@tiMainString  = title
forrest@0
   321
forrest@0
   322
  plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg) 
forrest@0
   323
forrest@0
   324
;(c) model-ob
forrest@0
   325
forrest@0
   326
  zz = datamod
forrest@0
   327
  zz = datamod - dataob
forrest@1
   328
  title = "Model "+model_name+" (2000) - Observed"
forrest@0
   329
forrest@0
   330
  resg@cnMinLevelValF  = -10.          ; Min level
forrest@0
   331
  resg@cnMaxLevelValF  =  10.          ; Max level
forrest@0
   332
  resg@cnLevelSpacingF =  2.           ; interval
forrest@0
   333
  resg@tiMainString    = title
forrest@0
   334
forrest@0
   335
  plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 
forrest@0
   336
forrest@0
   337
  pres                            = True        ; panel plot mods desired
forrest@0
   338
; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around
forrest@0
   339
                                                ; indiv. plots in panel
forrest@0
   340
  pres@gsnMaximize                = True        ; fill the page
forrest@0
   341
forrest@0
   342
  gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot
forrest@0
   343
forrest@0
   344
  delete (wks)
forrest@0
   345
  delete (plot)
forrest@0
   346
  delete (zz)
forrest@0
   347
forrest@0
   348
  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   349
         "rm "+plot_name+"."+plot_type)
forrest@1
   350
  ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0
   351
forrest@0
   352
  resg@gsnFrame             = True            ; draw plot 
forrest@0
   353
  resg@gsnDraw              = True            ; advance frame
forrest@0
   354
;------------------------------------------------------------------------
forrest@0
   355
; contour ob : masked
forrest@0
   356
 
forrest@0
   357
  resg@cnMinLevelValF       = 0.              ; Min level
forrest@0
   358
  resg@cnMaxLevelValF       = 30.             ; Max level
forrest@0
   359
  resg@cnLevelSpacingF      = 2.              ; interval
forrest@0
   360
  
forrest@0
   361
  plot_name = "global_mask_ob"
forrest@0
   362
  title     = ob_name+" "+ model_grid
forrest@0
   363
  resg@tiMainString  = title
forrest@0
   364
forrest@0
   365
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
   366
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
   367
;-----------------------------------------
forrest@0
   368
; plot average over mask region
forrest@0
   369
forrest@0
   370
  gRes  = True
forrest@0
   371
  gRes@txFontHeightF = 0.02
forrest@0
   372
  gRes@txAngleF = 0
forrest@0
   373
forrest@1
   374
  area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" kg C m~S~-2~N~)"
forrest@0
   375
forrest@0
   376
  gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
forrest@0
   377
;-----------------------------------------
forrest@0
   378
  zo = dataob
forrest@0
   379
  zo = dataob*mask_amazon
forrest@0
   380
  zo = where((mask_amazon .le. 0.01), zo@_FillValue, zo)  
forrest@0
   381
  plot = gsn_csm_contour_map_ce(wks,zo,resg)   
forrest@0
   382
forrest@0
   383
  delete (wks)
forrest@0
   384
  delete (plot)
forrest@0
   385
forrest@0
   386
  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   387
         "rm "+plot_name+"."+plot_type)
forrest@1
   388
  ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0
   389
forrest@0
   390
;------------------------------------------------------------------------
forrest@0
   391
; contour model: masked
forrest@0
   392
forrest@0
   393
  resg@cnMinLevelValF       = 0.              ; Min level
forrest@0
   394
  resg@cnMaxLevelValF       = 30.             ; Max level
forrest@0
   395
  resg@cnLevelSpacingF      = 2.              ; interval
forrest@0
   396
forrest@0
   397
  plot_name = "global_mask_model"
forrest@1
   398
  title     = "Model "+ model_name + " (2000)"
forrest@0
   399
  resg@tiMainString  = title
forrest@0
   400
forrest@0
   401
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
   402
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
   403
;-----------------------------------------
forrest@0
   404
; plot average over mask region
forrest@0
   405
forrest@0
   406
  gRes  = True
forrest@0
   407
  gRes@txFontHeightF = 0.02
forrest@0
   408
  gRes@txAngleF = 0
forrest@0
   409
forrest@1
   410
  area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" kg C m~S~-2~N~)"
forrest@0
   411
forrest@0
   412
  gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
forrest@0
   413
;-----------------------------------------
forrest@0
   414
  zm = datamod
forrest@0
   415
  zm = datamod*mask_amazon
forrest@0
   416
  zm = where((mask_amazon .le. 0.01), zm@_FillValue, zm)  
forrest@0
   417
  plot = gsn_csm_contour_map_ce(wks,zm,resg)     
forrest@0
   418
forrest@0
   419
  delete (wks)
forrest@0
   420
  delete (plot)
forrest@0
   421
forrest@0
   422
  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   423
         "rm "+plot_name+"."+plot_type)
forrest@1
   424
  ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0
   425
forrest@0
   426
;------------------------------------------------------------------------
forrest@0
   427
; contour model vs ob: masked 
forrest@0
   428
forrest@0
   429
  plot_name = "global_mask_vs_ob"
forrest@0
   430
forrest@0
   431
  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0
   432
  gsn_define_colormap(wks,"gui_default")     ; choose colormap
forrest@0
   433
forrest@0
   434
  plot=new(3,graphic)                        ; create graphic array
forrest@0
   435
forrest@0
   436
  resg@gsnFrame             = False          ; Do not draw plot 
forrest@0
   437
  resg@gsnDraw              = False          ; Do not advance frame
forrest@0
   438
forrest@0
   439
;(d) compute correlation coef and M score
forrest@0
   440
forrest@0
   441
  delete (good)
forrest@0
   442
  delete (uu)
forrest@0
   443
  delete (vv)
forrest@0
   444
  delete (ug)
forrest@0
   445
  delete (vg)  
forrest@0
   446
forrest@1
   447
  score_max = 0.
forrest@0
   448
forrest@0
   449
  uu = ndtooned(zm)
forrest@0
   450
  vv = ndtooned(zo)
forrest@0
   451
forrest@0
   452
  good = ind((uu .gt. 0.) .and. (vv .gt. 0.))
forrest@0
   453
forrest@0
   454
  ug = uu(good)
forrest@0
   455
  vg = vv(good)
forrest@0
   456
forrest@0
   457
  ccrG = esccr(ug,vg,0)
forrest@0
   458
forrest@0
   459
; Miomass = (ccrG*ccrG)*score_max 
forrest@0
   460
; new eq
forrest@0
   461
  bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
forrest@0
   462
  Mbiomass2  = (1. - (bias/dimsizes(ug)))*score_max
forrest@0
   463
  M_biomask = sprintf("%.2f", Mbiomass2)
forrest@0
   464
forrest@0
   465
  if (isvar("compare")) then
forrest@0
   466
     system("sed -e '1,/M_biomask/s/M_biomask/"+M_biomask+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0
   467
            "mv -f "+html_new2+" "+html_name2)
forrest@0
   468
     system("sed -e '1,/Sum_biomass_ob/s/Sum_biomass_ob/"+Sum_biomass_ob+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0
   469
            "mv -f "+html_new2+" "+html_name2)
forrest@0
   470
     system("sed -e '1,/Sum_biomass_mod/s/Sum_biomass_mod/"+Sum_biomass_mod+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0
   471
            "mv -f "+html_new2+" "+html_name2)
forrest@0
   472
  end if
forrest@0
   473
forrest@0
   474
  system("sed s#M_biomask#"+M_biomask+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   475
         "mv -f "+html_new+" "+html_name)
forrest@0
   476
  system("sed s#Sum_biomass_ob#"+Sum_biomass_ob+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   477
         "mv -f "+html_new+" "+html_name)
forrest@0
   478
  system("sed s#Sum_biomass_mod#"+Sum_biomass_mod+"# "+html_name+" > "+html_new+";"+ \
forrest@0
   479
         "mv -f "+html_new+" "+html_name)
forrest@0
   480
;--------------------------------------------------------------------
forrest@0
   481
; plot correlation coef
forrest@0
   482
forrest@0
   483
  gRes  = True
forrest@0
   484
  gRes@txFontHeightF = 0.02
forrest@0
   485
  gRes@txAngleF = 90
forrest@0
   486
forrest@0
   487
  correlation_text = "(correlation coef = "+sprintf("%.2f", ccrG)+")"
forrest@0
   488
forrest@0
   489
  gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
forrest@0
   490
;--------------------------------------------------------------------
forrest@0
   491
  
forrest@0
   492
;(a) ob
forrest@0
   493
forrest@0
   494
  title     = ob_name+" "+ model_grid
forrest@0
   495
  resg@tiMainString  = title
forrest@0
   496
forrest@0
   497
  plot(0) = gsn_csm_contour_map_ce(wks,zo,resg)       
forrest@0
   498
forrest@0
   499
;(b) model
forrest@0
   500
forrest@1
   501
  title     = "Model "+ model_name + " (2000)"
forrest@0
   502
  resg@tiMainString  = title
forrest@0
   503
forrest@0
   504
  plot(1) = gsn_csm_contour_map_ce(wks,zm,resg) 
forrest@0
   505
forrest@0
   506
;(c) model-ob
forrest@0
   507
forrest@0
   508
  zz = zo
forrest@0
   509
  zz = zm - zo
forrest@1
   510
  title = "Model "+model_name+" (2000) - Observed"
forrest@0
   511
forrest@0
   512
  resg@cnMinLevelValF  = -10.          ; Min level
forrest@0
   513
  resg@cnMaxLevelValF  =  10.          ; Max level
forrest@0
   514
  resg@cnLevelSpacingF =  2.           ; interval
forrest@0
   515
  resg@tiMainString    = title
forrest@0
   516
forrest@0
   517
  plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 
forrest@0
   518
forrest@0
   519
  pres                            = True        ; panel plot mods desired
forrest@0
   520
; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around
forrest@0
   521
                                                ; indiv. plots in panel
forrest@0
   522
  pres@gsnMaximize                = True        ; fill the page
forrest@0
   523
forrest@0
   524
  gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot
forrest@0
   525
forrest@0
   526
  delete (wks)
forrest@0
   527
  delete (plot)
forrest@0
   528
forrest@0
   529
  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
forrest@0
   530
         "rm "+plot_name+"."+plot_type)
forrest@1
   531
  ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0
   532
forrest@0
   533
;***************************************************************************
forrest@0
   534
; add total score and write to file
forrest@0
   535
;***************************************************************************
forrest@0
   536
  M_total = Mbiomass + Mbiomass2
forrest@0
   537
forrest@0
   538
  asciiwrite("M_save.biomass", M_total)
forrest@0
   539
forrest@0
   540
;***************************************************************************
forrest@0
   541
; output plots
forrest@0
   542
;***************************************************************************
forrest@0
   543
  output_dir = model_name+"/biomass"
forrest@0
   544
forrest@0
   545
  system("mv *.png " + output_dir) 
forrest@0
   546
;***************************************************************************
forrest@0
   547
end