var/graphics/ncl/gen_be/gsi_be_plots.ncl

   1 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
   4
   5    ; just edit the filename for your case
   6
   7    filename = "/Volumes/sugar2/hclin/data/tmp_afwa/gen_be/run_gsi_be_lat_bin_size_5.0_lnps/wrf-arw-gsi_be.gcv"
   8    ;filename = "/Volumes/sugar2/hclin/code/GSI/trunk/fix/nam_nmm_berror.f77.gcv"
   9    ;filename = "/Volumes/sugar2/hclin/code/GSI/trunk/fix/nam_glb_berror.f77.gcv"
  10
  11    plot_dir = "."
  12    out_name = "gsi_be_plots"
  13    out_type = "pdf"
  14
  15    Fill_Color = True
  16    By_Levels = True
  17    sig_values_for_plot=(/0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1/) ; settings for By_Levels = False
  18    sig_array_for_plot=(/" .9"," .8"," .7"," .6"," .5"," .4"," .3"," .2"," .1"/) ; settings for By_Levels = False
  19    ;tmMode = "Explicit"
  20    tmMode = "Automatic"
  21    lat_values_for_plot =(/ 0, 10, 20, 30, 40, 50, 60, 70/)  ; settings for tmMode = "Explicit"
  22    lat_array_for_plot = (/"0N","10N","20N","30N","40N","50N","60N","70N"/)  ; ; settings for tmMode = "Explicit"
  23    plot_dims_aero = (/5,4/)
  24
  25    varnames_met = (/ "sf", "vp", "t", "q", "ps" /)
  26    nvar_met = dimsizes(varnames_met)
  27
  28    varnames_aero = (/ "BC1", "BC2", "OC1", "OC2", "SEAS_1", "SEAS_2", "SEAS_3", "SEAS_4", \
  29                       "DUST_1", "DUST_2", "DUST_3", "DUST_4", "DUST_5", "sulf", "P25" /)
  30    nvar_aero = dimsizes(varnames_aero)
  31
  32    if ( .not. isfilepresent(filename) ) then
  33       print("Error: can not find "+filename)
  34       exit
  35    end if
  36
  37    ;-----------------
  38    ; read in data
  39    ;-----------------
  40    setfileoption("bin","ReadByteOrder","BigEndian")
  41    dims = fbinrecread(filename,0,2,"integer")
  42    nsig=dims(0)
  43    nlat=dims(1)
  44    print("nsig, nlat = "+nsig+", "+nlat)
  45
  46    bytes_in_file = stringtointeger(systemfunc("wc -c "+filename))
  47
  48    ; one record = 4 bytes
  49    ; one character = 1 byte
  50    ; end of record = 8 bytes
  51    nrec_met = 3+4*(nvar_met-1)+3
  52    ncount_met = 2+(nlat+nsig)+(nsig*nsig*(nlat+2))+(nsig*(nlat+2))+(nsig*(nlat+2))+ \
  53                 (nlat*nsig)*nvar_met+nlat+ \
  54                 nlat+2+(nvar_met-1)*((nlat+2)*nsig+(nlat+2)*nsig) +nvar_met
  55    nbyte_met = ncount_met*4+5*nvar_met+8*nrec_met
  56
  57    nrec_aero = 4*nvar_aero
  58    ncount_aero = nvar_aero*(1+(nlat*nsig)+(nlat+2)*nsig+(nlat+2)*nsig)
  59    nbyte_aero = ncount_aero*4+5*nvar_aero+8*nrec_aero
  60    nbyte_tot = nbyte_met+nbyte_aero
  61
  62    read_aero = False
  63    nrec_tot = nrec_met
  64    nvar_tot = nvar_met
  65    varnames_tot = varnames_met
  66    if ( bytes_in_file .eq. nbyte_tot ) then
  67       read_aero = True
  68       nrec_tot = nrec_met + nrec_aero
  69       nvar_tot = nvar_met + nvar_aero
  70       delete(varnames_tot)
  71       varnames_tot = array_append_record (varnames_met, varnames_aero, 0)
  72    end if
  73
  74    rlat = new(nlat, "float")
  75    rsig = new(nsig, "float")
  76    agv  = new((/nsig,nsig,nlat+2/), "float")
  77    bv   = new((/nsig,nlat+2/), "float")
  78    wgv  = new((/nsig,nlat+2/), "float")
  79    cov  = new((/nvar_tot,nsig,nlat/), "float")
  80    covq2 = new((/nsig,nlat/), "float")
  81    ;hwllp = new(nlat+2, "float")
  82    hwll  = new((/nvar_tot,nsig,nlat+2/), "float")
  83    ;vzs   = new((/nvar_tot,nsig,nlat+2/), "float")
  84    vzs   = new((/nvar_tot,nlat+2,nsig/), "float")
  85
  86    rtmp = fbinrecread(filename,1,(nlat+nsig),"float")
  87    rlat = rtmp(0:nlat-1)
  88    rsig = rtmp(nlat:nlat+nsig-1)
  89    delete(rtmp)
  90
  91    yc = ispan(1,nsig,1)
  92    rlatp2 = new(nlat+2, "float")
  93    rlatp2(0) = rlat(0)
  94    rlatp2(nlat+1) = rlat(nlat-1)
  95    rlatp2(1:nlat) = rlat(0:nlat-1)
  96    ; do k = 0, nlat-2
  97    ;   rlatp2(k+1) = 0.5*(rlat(k)+rlat(k+1))
  98    ;end do
  99    ;rlatp2(nlat)   = 0.5*(rlatp2(nlat-1)+rlat(nlat-1))
 100    ;rlatp2(nlat+1) = rlat(nlat-1)
 101
 102    dsig = new(nsig, "float")
 103    dsig(0)=log(rsig(0))-log(rsig(1))
 104    do k=1,nsig-2
 105       dsig(k)=0.5*(log(rsig(k-1))-log(rsig(k+1)))
 106    end do
 107    dsig(nsig-1)=log(rsig(nsig-2))-log(rsig(nsig-1))
 108
 109    cov!0 = "var"
 110    cov!1 = "lev"
 111    cov!2 = "lat"
 112    cov&var = varnames_tot
 113    cov&lat = rlat
 114
 115    hwll!0 = "var"
 116    hwll!1 = "lev"
 117    hwll!2 = "lat"
 118    hwll&var = varnames_tot
 119    hwll&lat = rlatp2
 120
 121    vzs!0 = "var"
 122    vzs!1 = "lat"
 123    vzs!2 = "lev"
 124    vzs&var = varnames_tot
 125    vzs&lat = rlatp2
 126
 127    if ( .not. By_Levels ) then
 128       var&lev = rsig
 129       hwll&lev = rsig
 130       vzs&lev = rsig
 131    else
 132       cov&lev = yc
 133       hwll&lev = yc
 134       vzs&lev = yc
 135    end if
 136
 137    reclen = (nsig*nsig*(nlat+2))+(nsig*(nlat+2))+(nsig*(nlat+2))
 138    rtmp = fbinrecread(filename,2,reclen,"float")
 139    is = 0
 140    ie = nsig*nsig*(nlat+2)-1
 141    agv = onedtond(rtmp(is:ie),(/nsig,nsig,nlat+2/))
 142    is = nsig*nsig*(nlat+2)
 143    ie = (nsig*nsig*(nlat+2))+(nsig*(nlat+2))-1
 144    bv  = onedtond(rtmp(is:ie),(/nsig,nlat+2/))
 145    is = (nsig*nsig*(nlat+2))+(nsig*(nlat+2))
 146    ie = (nsig*nsig*(nlat+2))+(nsig*(nlat+2))+(nsig*(nlat+2))-1
 147    wgv  = onedtond(rtmp(is:ie),(/nsig,nlat+2/))
 148    delete(rtmp)
 149
 150    do i = 1, nvar_met
 151       irec = 4*(i-1)+3
 152       rtmp = fbinrecread(filename,irec,9,"character")  ; 5 bytes plus 4 bytes
 153       var = chartostring(rtmp(0:4))
 154       print("Reading "+var)
 155       delete(rtmp)
 156
 157       irec = 4*(i-1)+4
 158       if ( var .eq. "q    " ) then
 159          reclen = 2*nlat*nsig
 160       else
 161          if ( var .ne. "ps   " ) then
 162             reclen = nlat*nsig
 163          else
 164             reclen = nlat
 165          end if
 166       end if
 167       rtmp = fbinrecread(filename,irec,reclen,"float")
 168       if ( var .eq. "q    " ) then
 169          cov(i-1,:,:) = onedtond(rtmp(0:nlat*nsig-1),(/nsig,nlat/))
 170          covq2(:,:) = onedtond(rtmp(nlat*nsig:2*nlat*nsig-1),(/nsig,nlat/))
 171          cov(i-1,:,:) = cov(i-1,:,:)*100.0
 172          covq2(:,:) = covq2(:,:)*100.0
 173       else
 174          if ( var .ne. "ps   " ) then
 175             cov(i-1,:,:) = onedtond(rtmp,(/nsig,nlat/))
 176             if ( var .eq. "sf   " .or. \
 177                  var .eq. "vp   " ) then
 178                cov(i-1,:,:) = cov(i-1,:,:) * 0.000001
 179             end if
 180          else
 181             cov(i-1,0,:) = (/rtmp/)
 182          end if
 183       end if
 184       delete(rtmp)
 185
 186       if ( var .ne. "ps   " ) then
 187          irec = 4*(i-1)+5
 188          reclen = (nlat+2)*nsig
 189          rtmp = fbinrecread(filename,irec,reclen,"float")
 190          hwll(i-1,:,:) = onedtond(rtmp,(/nsig,nlat+2/))
 191          delete(rtmp)
 192          irec = 4*(i-1)+6
 193          reclen = (nlat+2)*nsig
 194          rtmp = fbinrecread(filename,irec,reclen,"float")
 195          vzs(i-1,:,:) = onedtond(rtmp,(/nlat+2,nsig/))
 196          delete(rtmp)
 197       else
 198          irec = 4*(i-1)+5
 199          reclen = nlat+2
 200          rtmp = fbinrecread(filename,irec,reclen,"float")
 201          hwll(i-1,0,:) = (/rtmp/)
 202          delete(rtmp)
 203       end if
 204    end do
 205
 206    if ( read_aero ) then
 207       do i = nvar_met+1, nvar_tot
 208          irec = 4*(i-1)+2
 209          rtmp = fbinrecread(filename,irec,9,"character")  ; 5 bytes plus 4 bytes
 210          var = chartostring(rtmp(0:4))
 211          print("Reading "+var)
 212          delete(rtmp)
 213
 214          irec = 4*(i-1)+3
 215          reclen = nlat*nsig
 216          rtmp = fbinrecread(filename,irec,reclen,"float")
 217          cov(i-1,:,:) = onedtond(rtmp,(/nsig,nlat/))
 218          delete(rtmp)
 219
 220          irec = 4*(i-1)+4
 221          reclen = (nlat+2)*nsig
 222          rtmp = fbinrecread(filename,irec,reclen,"float")
 223          hwll(i-1,:,:) = onedtond(rtmp,(/nsig,nlat+2/))
 224          delete(rtmp)
 225          irec = 4*(i-1)+5
 226          reclen = (nlat+2)*nsig
 227          rtmp = fbinrecread(filename,irec,reclen,"float")
 228          vzs(i-1,:,:) = onedtond(rtmp,(/nlat+2,nsig/))
 229          delete(rtmp)
 230       end do
 231    end if
 232
 233    hwll    = hwll*0.001 ; km
 234    do i = 0, nvar_tot-1
 235       do j = 0, nlat+1
 236          do k = 0, nsig-1
 237             if ( .not. ismissing(vzs(i,j,k)) ) then
 238                if ( vzs(i,j,k) .gt. 0.0 ) then
 239                   vzs(i,j,k) = 1.0/vzs(i,j,k)/dsig(k)
 240                else
 241                   vzs(i,j,k) = 0.0
 242                end if
 243             end if
 244          end do
 245       end do
 246    end do
 247    reg_t   = agv*1000000000.0
 248    reg_chi = bv
 249    reg_ps  = wgv*1000000000.0
 250
 251    reg_t!0   = "lev"
 252    reg_t!1   = "lev"
 253    reg_t!2   = "lat"
 254    reg_t&lat = rlatp2
 255
 256    reg_chi!0   = "lev"
 257    reg_chi!1   = "lat"
 258    reg_chi&lat = rlatp2
 259
 260    reg_ps!0   = "lev"
 261    reg_ps!1   = "lat"
 262    reg_ps&lat = rlatp2
 263
 264    if ( .not. By_Levels ) then
 265       reg_t&lev   = rsig
 266       reg_ps&lev  = rsig
 267       reg_chi&lev = rsig
 268    else
 269       reg_t&lev   = yc
 270       reg_ps&lev  = yc
 271       reg_chi&lev = yc
 272    end if
 273
 274    ; averages (over latitudes) for each level
 275    cov_avg  = dim_avg_n_Wrap(cov, 2)  ; (nvar, nsig, nlat)->(nvar,nsig)
 276    hwll_avg = dim_avg_n_Wrap(hwll, 2) ; (nvar, nsig, nlat+2)->(nvar,nsig)
 277    vzs_avg  = dim_avg_n_Wrap(vzs, 1)  ; (nvar, nlat+2, nsig)->(nvar,nsig)
 278
 279    ;----------------------
 280    ; start plotting
 281    ;----------------------
 282    wks = gsn_open_wks (out_type,plot_dir+"/"+out_name)  ; open workstation
 283    gsn_define_colormap(wks,"rainbow+white+gray")
 284
 285    ;--------------------
 286    ; Plot Variance
 287    ;--------------------
 288    print(" 1. Plotting Standard Deviation")
 289    plts = new (nvar_tot,"graphic")
 290
 291    res                       = True
 292    res@cnFillOn              = Fill_Color
 293    res@gsnSpreadColors       = True
 294    res@gsnSpreadColorStart   = 20
 295    res@gsnSpreadColorEnd     = -2
 296    res@lbLabelAutoStride     = True
 297    res@lbOrientation         = "Vertical"
 298    res@gsnDraw               = False
 299    res@gsnFrame              = False
 300    res@tmXTOn                = False
 301    res@tmYROn                = False
 302    if ( .not. By_Levels ) then
 303       res@tiYAxisString = "Sigma Values"
 304       res@trYReverse = True
 305       res@tmYLMode   = "Explicit"
 306       res@tmYLValues = sig_values_for_plot
 307       res@tmYLLabels = sig_array_for_plot
 308    else
 309       res@tiYAxisString = "Sigma Levels"
 310    end if
 311    res@tiXAxisString = "Latitude"
 312    res@tmXBMode         = tmMode
 313    if( tmMode .eq. "Explicit") then
 314       res@tmXBValues    = lat_values_for_plot
 315       res@tmXBLabels    = lat_array_for_plot
 316    end if
 317    res@gsnYAxisIrregular2Linear=True
 318
 319    resP                     = True
 320    resP@gsnMaximize         = True
 321    resP@gsnPaperOrientation = "portrait"
 322    resP@txString            = "Standard Deviation"
 323    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 324
 325    do kk = 0,nvar_met-2
 326       res@gsnLeftString = str_squeeze(varnames_met(kk))
 327       res@gsnRightString = ""
 328       if ( kk .eq. 0 .or. kk .eq. 1 ) then  ; sf, vp
 329          res@gsnRightString = "x 10**6"
 330       end if
 331       if ( kk .eq. 3 ) then  ; q
 332          res@gsnRightString = "x 0.01"
 333       end if
 334       plts(kk) = gsn_csm_contour (wks,cov(kk,:,:),res)
 335    end do
 336    gsn_panel(wks,plts(0:nvar_met-2),(/2,2/),resP)
 337
 338    if ( read_aero ) then
 339       delete(res@gsnRightString)
 340       do kk = 0,nvar_aero-1
 341          res@gsnLeftString = str_squeeze(varnames_aero(kk))
 342          plts(kk+nvar_met) = gsn_csm_contour (wks,cov(kk+nvar_met,:,:),res)
 343       end do
 344       gsn_panel(wks,plts(nvar_met:nvar_tot-1),plot_dims_aero,resP)
 345    end if
 346
 347    delete (res)
 348    delete (resP)
 349    delete (plts)
 350
 351    ;--------------------------------
 352    ; Plot Average Variance Profile
 353    ;--------------------------------
 354    print(" 2. Plotting Averaged Standard Deviation Profile")
 355    plts = new (nvar_tot,"graphic")
 356
 357    res                    = True
 358    res@gsnDraw            = False
 359    res@gsnFrame           = False
 360    res@tiYAxisString      = "Sigma Values"
 361    res@xyLineThicknesses  = (/4.0,4.0,4.0,4.0/)
 362    res@xyDashPatterns     = (/0,4,2,3/)
 363    res@xyLineColors       = (/"blue"/)
 364    res@tmXTOn             = False
 365    res@tmYROn             = False
 366    if ( .not. By_Levels ) then
 367       res@tiYAxisString = "Sigma Values"
 368       res@trYReverse    = True
 369       res@tmYLMode      = "Explicit"
 370       res@tmYLValues    = sig_values_for_plot
 371       res@tmYLLabels    = sig_array_for_plot
 372    else
 373       res@trYMinF  =  1.0
 374       res@trYMaxF  =  nsig + 0.0
 375       res@tiYAxisString   = "Sigma Levels"
 376    end if
 377
 378    res@lgPerimOn              = False
 379    res@lgLabelFontHeightF     = 0.02
 380    res@pmLegendDisplayMode    = "Always"
 381    res@pmLegendSide           = "Bottom"
 382    res@pmLegendParallelPosF   = 0.20
 383    res@pmLegendOrthogonalPosF = -0.45
 384    res@pmLegendWidthF         = 0.2
 385    res@pmLegendHeightF        = 0.2
 386
 387    resP                     = True
 388    resP@gsnMaximize         = True
 389    resP@gsnPaperOrientation = "portrait"
 390    resP@txString            = "Averaged Standard Deviation Profile"
 391    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 392
 393    do kk = 0, nvar_met-2
 394       res@xyExplicitLegendLabels = varnames_met(kk)
 395       res@tiXAxisString   = "Standard Deviation"
 396       if ( .not. By_Levels ) then
 397          plts(kk) = gsn_csm_xy (wks,cov_avg(kk,:),rsig,res)
 398       else
 399          plts(kk) = gsn_csm_xy (wks,cov_avg(kk,:),yc,res)
 400       end if
 401    end do
 402    gsn_panel(wks,plts(0:nvar_met-2),(/2,2/),resP)
 403
 404    if ( read_aero ) then
 405       do kk = 0, nvar_aero-1
 406          res@xyExplicitLegendLabels = varnames_aero(kk)
 407          res@tiXAxisString   = "Standard Deviation"
 408          if ( .not. By_Levels ) then
 409             plts(kk+nvar_met) = gsn_csm_xy (wks,cov_avg(kk+nvar_met,:),rsig,res)
 410          else
 411             plts(kk+nvar_met) = gsn_csm_xy (wks,cov_avg(kk+nvar_met,:),yc,res)
 412          end if
 413       end do
 414       gsn_panel(wks,plts(nvar_met:nvar_tot-1),plot_dims_aero,resP)
 415    end if
 416
 417    delete (res )
 418    delete (resP )
 419    delete (plts)
 420
 421    ;-----------------------------
 422    ; Plot Horizontal len scales
 423    ;-----------------------------
 424    print(" 3. Plotting Horizontal Length-scale")
 425    plts = new (nvar_tot,"graphic")
 426
 427    res                       = True
 428    res@cnFillOn              = Fill_Color
 429    res@gsnSpreadColors       = True
 430    res@gsnSpreadColorStart   = 20
 431    res@gsnSpreadColorEnd     = -2
 432    res@lbLabelAutoStride     = True
 433    res@lbOrientation         = "Vertical"
 434    res@gsnDraw               = False
 435    res@gsnFrame              = False
 436    res@tmXTOn                = False
 437    res@tmYROn                = False
 438    res@gsnYAxisIrregular2Linear=True
 439    if ( .not. By_Levels ) then
 440       res@tiYAxisString = "Sigma Values"
 441       res@trYReverse = True
 442       res@tmYLMode   = "Explicit"
 443       res@tmYLValues = sig_values_for_plot
 444       res@tmYLLabels = sig_array_for_plot
 445    else
 446       res@tiYAxisString = "Sigma Levels"
 447    end if
 448    res@tiXAxisString = "Latitude"
 449    res@tmXBMode = tmMode
 450    if ( res@tmXBMode .eq. "Explicit" ) then
 451       res@tmXBValues = lat_values_for_plot
 452       res@tmXBLabels = lat_array_for_plot
 453    end if
 454
 455    resP                     = True
 456    resP@gsnMaximize         = True
 457    resP@gsnPaperOrientation = "portrait"
 458    resP@txString            = "Horizontal Length-scale (Km)"
 459    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 460
 461    do kk = 0, nvar_met-2
 462       res@gsnLeftString = str_squeeze(varnames_met(kk))
 463       plts(kk) = gsn_csm_contour (wks,hwll(kk,:,1:nlat),res)
 464    end do
 465    gsn_panel(wks,plts(0:nvar_met-2),(/2,2/),resP)
 466
 467    if ( read_aero ) then
 468       do kk = 0, nvar_aero-1
 469          res@gsnLeftString = str_squeeze(varnames_aero(kk))
 470          plts(kk+nvar_met) = gsn_csm_contour (wks,hwll(kk+nvar_met,:,1:nlat),res)
 471       end do
 472       gsn_panel(wks,plts(nvar_met:nvar_tot-1),plot_dims_aero,resP)
 473    end if
 474
 475    delete (res)
 476    delete (resP)
 477    delete (plts)
 478
 479    ;-------------------------------------
 480    ; Plot Avergaed Horizontal len scales
 481    ;-------------------------------------
 482    print(" 4. Plotting Averaged Horizontal Length-scale Profile")
 483    plts = new (nvar_tot,"graphic")
 484
 485    res                    = True
 486    res@gsnDraw            = False
 487    res@gsnFrame           = False
 488    res@tiXAxisString      = "Horizontal Length-scale (Km)"
 489    res@xyLineThicknesses  = (/4.0,4.0,4.0,4.0/)
 490    res@xyLineColors       = (/"blue","blue","green","purple"/)
 491    res@xyDashPatterns     = (/0,4,2,3,4/)
 492    res@tmXTOn             = False
 493    res@tmYROn             = False
 494    if ( .not. By_Levels ) then
 495       res@tiYAxisString = "Sigma Values"
 496       res@trYReverse    = True
 497       res@tmYLMode      = "Explicit"
 498       res@tmYLValues    = sig_values_for_plot
 499       res@tmYLLabels    = sig_array_for_plot
 500    else
 501       res@trYMinF  =  1.0
 502       res@trYMaxF  =  nsig + 0.0
 503       res@tiYAxisString   = "Sigma Levels"
 504    end if
 505
 506    res@lgPerimOn              = False
 507    res@lgLabelFontHeightF     = 0.02
 508    res@pmLegendDisplayMode    = "Always"
 509    res@pmLegendSide           = "Bottom"
 510    res@pmLegendParallelPosF   = 0.83
 511    res@pmLegendOrthogonalPosF = -0.40
 512    res@pmLegendWidthF         = 0.2
 513    res@pmLegendHeightF        = 0.2
 514
 515    resP                     = True
 516    resP@gsnMaximize         = True
 517    resP@gsnPaperOrientation = "portrait"
 518    resP@txString            = "Averaged Horizontal Length-scale (km)"
 519    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 520
 521    do kk = 0, nvar_met-2
 522       res@xyExplicitLegendLabels = varnames_met(kk)
 523       if ( .not. By_Levels ) then
 524          plts(kk) = gsn_csm_xy (wks,hwll_avg(kk,:),rsig,res)
 525       else
 526          plts(kk) = gsn_csm_xy (wks,hwll_avg(kk,:),yc,res)
 527       end if
 528    end do
 529    gsn_panel(wks,plts(0:nvar_met-2),(/2,2/),resP)
 530
 531    if ( read_aero ) then
 532       do kk = 0, nvar_aero-1
 533          res@xyExplicitLegendLabels = varnames_aero(kk)
 534          if ( .not. By_Levels ) then
 535             plts(kk+nvar_met) = gsn_csm_xy (wks,hwll_avg(kk+nvar_met,:),rsig,res)
 536          else
 537             plts(kk+nvar_met) = gsn_csm_xy (wks,hwll_avg(kk+nvar_met,:),yc,res)
 538          end if
 539       end do
 540       gsn_panel(wks,plts(nvar_met:nvar_tot-1),plot_dims_aero,resP)
 541    end if
 542
 543    delete (res)
 544    delete (resP)
 545    delete (plts)
 546
 547    ;---------------
 548    ; Plot Ps
 549    ;---------------
 550    print(" 5. Plotting Variance & Horizontal Length-scale for surface pressure")
 551    plts = new (2,"graphic")
 552
 553    res                   = True
 554    res@gsnDraw           = False
 555    res@gsnFrame          = False
 556    res@xyLineThicknesses = (/4.0,4.0,4.0,4.0/)
 557    res@xyDashPatterns    = (/0,4,2,3/)
 558    res@xyLineColors      = (/"blue"/)
 559    res@tmXTOn            = False
 560    res@tmYROn            = False
 561    res@tiXAxisString   = "Latitude"
 562    res@tmXBMode = tmMode
 563    if ( res@tmXBMode .eq. "Explicit" ) then
 564       res@tmXBValues = lat_values_for_plot
 565       res@tmXBLabels = lat_array_for_plot
 566    end if
 567    res@tiYAxisString   = "Standard Deviation (hPa)"
 568
 569    res@lgPerimOn              = False
 570    res@lgLabelFontHeightF     = 0.02
 571    res@pmLegendDisplayMode    = "Always"
 572    res@pmLegendSide           = "Bottom"
 573    res@pmLegendParallelPosF   = 0.75
 574    res@pmLegendOrthogonalPosF = -0.45
 575    res@pmLegendWidthF         = 0.2
 576    res@pmLegendHeightF        = 0.2
 577    res@xyExplicitLegendLabels = "ps_u"
 578
 579    plts(0) = gsn_csm_xy(wks,rlat,cov(4,0,:),res)
 580
 581    res@tiYAxisString   = "Lengthscale (Km)"
 582    plts(1) = gsn_csm_xy(wks,rlat,hwll(4,0,1:nlat),res)
 583
 584    resP                     = True
 585    resP@gsnMaximize         = True
 586    resP@gsnPaperOrientation = "portrait"
 587    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 588    gsn_panel(wks,plts,(/2,1/),resP)
 589
 590    delete (res)
 591    delete (resP)
 592    delete (plts)
 593
 594    ;--------------------------
 595    ; Plot Vertical len scales
 596    ;--------------------------
 597    print(" 6. Plotting Vertical Length-scale ")
 598    plts  = new (nvar_tot,"graphic")
 599
 600    res                       = True
 601    res@gsnDraw               = False
 602    res@gsnFrame              = False
 603    res@cnFillOn              = Fill_Color
 604    res@gsnSpreadColors       = True
 605    res@gsnSpreadColorStart   = 20
 606    res@gsnSpreadColorEnd     = -2
 607    res@lbLabelAutoStride     = True
 608    res@lbOrientation         = "Vertical"
 609    res@tmXTOn                = False
 610    res@tmYROn                = False
 611    if ( .not. By_Levels ) then
 612       res@tiYAxisString   = "Sigma Values"
 613       res@trYReverse        = True
 614       res@tmYLMode = "Explicit"
 615       res@tmYLValues = sig_values_for_plot
 616       res@tmYLLabels = sig_array_for_plot
 617    else
 618       res@tiYAxisString   = "Sigma Levels"
 619    end if
 620    res@tiXAxisString = "Latitude"
 621    res@tmXBMode = tmMode
 622    if ( res@tmXBMode .eq. "Explicit" ) then
 623       res@tmXBMode      = "Explicit"
 624       res@tmXBValues    = lat_values_for_plot
 625       res@tmXBLabels    = lat_array_for_plot
 626    end if
 627
 628    resP                     = True
 629    resP@gsnMaximize         = True
 630    resP@gsnPaperOrientation = "portrait"
 631    resP@txString            = "Vertical Length-scale (sigma units)"
 632    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 633
 634    do kk = 0, nvar_met-2
 635       res@gsnLeftString         = str_squeeze(varnames_met(kk))
 636       plts(kk) = gsn_csm_contour(wks,vzs(var|kk,lev|:,lat|1:nlat),res)
 637    end do
 638    gsn_panel(wks,plts(0:nvar_met-2),(/2,2/),resP)
 639
 640    if ( read_aero ) then
 641       do kk = 0, nvar_aero-1
 642          res@gsnLeftString = str_squeeze(varnames_aero(kk))
 643          plts(kk+nvar_met) = gsn_csm_contour(wks,vzs(var|(kk+nvar_met),lev|:,lat|1:nlat),res)
 644       end do
 645       gsn_panel(wks,plts(nvar_met:nvar_tot-1),plot_dims_aero,resP)
 646    end if
 647
 648    delete (res)
 649    delete (resP)
 650    delete (plts)
 651
 652    ;-----------------------------------
 653    ; Plot Avergaed Vertical len scales
 654    ;-----------------------------------
 655    print(" 7. Plotting Averaged Vertical Length-scale Profile")
 656    plts = new (nvar_tot,"graphic")
 657
 658    res                   = True
 659    res@gsnDraw           = False
 660    res@gsnFrame          = False
 661    res@xyLineThicknesses = (/4.0,4.0,4.0,4.0/)
 662    res@xyLineColors      = (/"blue","blue","green","purple"/)
 663    res@xyDashPatterns    = (/0,4,2,3,4/)
 664    res@tmXTOn            = False
 665    res@tmYROn            = False
 666    if ( .not. By_Levels ) then
 667       res@tiYAxisString  = "Sigma Values"
 668       res@trYReverse = True
 669       res@tmYLMode   = "Explicit"
 670       res@tmYLValues = sig_values_for_plot
 671       res@tmYLLabels = sig_array_for_plot
 672    else
 673       res@tiYAxisString = "Sigma Levels"
 674       res@trYMinF       =  1.0
 675       res@trYMaxF       =  nsig + 0.0
 676    end if
 677    res@tiXAxisString     = "Vertical Lengthscale (sigma units)"
 678
 679    res@lgPerimOn              = False
 680    res@lgLabelFontHeightF     = 0.02
 681    res@pmLegendDisplayMode    = "Always"
 682    res@pmLegendSide           = "Bottom"
 683    res@pmLegendParallelPosF   = 0.85
 684    res@pmLegendOrthogonalPosF = -1.25
 685    res@pmLegendWidthF         = 0.2
 686    res@pmLegendHeightF        = 0.2
 687
 688    resP                     = True
 689    resP@gsnMaximize         = True
 690    resP@gsnPaperOrientation = "portrait"
 691    resP@txString            = "Averaged Vertical Length-scale (sigma units)"
 692    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 693
 694    do kk = 0, nvar_met-2
 695       res@xyExplicitLegendLabels = varnames_met(kk)
 696       if ( .not. By_Levels ) then
 697          plts(kk) = gsn_csm_xy (wks,vzs_avg(kk,:),rsig,res)
 698       else
 699          plts(kk) = gsn_csm_xy (wks,vzs_avg(kk,:),yc,res)
 700       end if
 701    end do
 702    gsn_panel(wks,plts(0:nvar_met-2),(/2,2/),resP)
 703
 704    if ( read_aero ) then
 705       do kk = 0, nvar_aero-1
 706          res@xyExplicitLegendLabels = varnames_aero(kk)
 707          if ( .not. By_Levels ) then
 708             plts(kk+nvar_met) = gsn_csm_xy (wks,vzs_avg(kk+nvar_met,:),rsig,res)
 709          else
 710             plts(kk+nvar_met) = gsn_csm_xy (wks,vzs_avg(kk+nvar_met,:),yc,res)
 711          end if
 712       end do
 713       gsn_panel(wks,plts(nvar_met:nvar_tot-1),plot_dims_aero,resP)
 714    end if
 715
 716    delete (res)
 717    delete (resP)
 718    delete (plts)
 719
 720    ;-----------------------------------
 721    ; Plot Chi Regression coefficients
 722    ;-----------------------------------
 723    print(" 8. Plotting Chi Regression coefficients")
 724    plts = new (2,"graphic")
 725
 726    res                       = True
 727    res@gsnMaximize           = True
 728    res@gsnDraw               = False
 729    res@gsnFrame              = False
 730    res@tmXTOn                = False
 731    res@tmYROn                = False
 732    res@cnFillOn              = Fill_Color
 733    res@lbLabelAutoStride     = True
 734    res@lbOrientation         = "Vertical"
 735    res@gsnSpreadColors       = True
 736    res@gsnSpreadColorStart   = 20
 737    res@gsnSpreadColorEnd     = -2
 738    res@gsnYAxisIrregular2Linear=True
 739    if ( .not. By_Levels ) then
 740       res@tiYAxisString = "Sigma Values"
 741       res@trYReverse = True
 742       res@tmYLMode   = "Explicit"
 743       res@tmYLValues = sig_values_for_plot
 744       res@tmYLLabels = sig_array_for_plot
 745    else
 746       res@tiYAxisString = "Sigma Levels"
 747       res@trYMinF       =  1.0
 748       res@trYMaxF       =  nsig + 0.0
 749    end if
 750    res@tiXAxisString = "Latitude"
 751    res@tmXBMode = tmMode
 752    if ( res@tmXBMode .eq. "Explicit" ) then
 753       res@tmXBValues = lat_values_for_plot
 754       res@tmXBLabels = lat_array_for_plot
 755    end if
 756    res@gsnLeftString   = "Chi - Regression Coefficients"
 757    plts(0) = gsn_csm_contour (wks,reg_chi(:,1:nlat),res)
 758    delete (res)
 759
 760    ;--------------------------------------------------------------------
 761    ;  Plot averaged Vertical profile of Chi regression coeff
 762    ;--------------------------------------------------------------------
 763
 764    reg_chi_avg = dim_avg_n_Wrap(reg_chi,1)
 765
 766    res                   = True
 767    res@gsnMaximize       = True
 768    res@gsnDraw           = False
 769    res@gsnFrame          = False
 770    res@tmXTOn            = False
 771    res@tmYROn            = False
 772    res@xyLineThicknesses = (/4.0,4.0/)
 773    res@xyDashPatterns    = (/0,4/)
 774    res@xyLineColors      = (/"blue"/)
 775
 776    res@tiXAxisString = "Regression Coefficient"
 777    if ( .not. By_Levels ) then
 778       res@tiYAxisString = "Sigma Values"
 779       res@trYReverse    = True
 780       res@tmYLMode      = "Explicit"
 781       res@tmYLValues    = sig_values_for_plot
 782       res@tmYLLabels    = sig_array_for_plot
 783    else
 784       res@tiYAxisString = "Sigma Levels"
 785       res@trYMinF       =  1.0
 786       res@trYMaxF       =  nsig + 0.0
 787    end if
 788
 789    res@xyExplicitLegendLabels = "chi_u"
 790    res@gsnLeftString   = "Chi - Averaged Regression Coefficients"
 791    if ( .not. By_Levels ) then
 792       plts(1) = gsn_csm_xy (wks,reg_chi_avg,rsig,res)
 793    else
 794       plts(1) = gsn_csm_xy (wks,reg_chi_avg,yc,res)
 795    end if
 796
 797    resP                     = True
 798    resP@gsnMaximize         = True
 799    resP@gsnPaperOrientation = "portrait"
 800    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 801    resP@gsnPanelBottom      = 0.10
 802    gsn_panel(wks,plts,(/2,1/),resP)
 803
 804    delete (res)
 805    delete (resP)
 806    delete (plts)
 807
 808    ;------------------------------------
 809    ; Plot Psfc Regression coefficients
 810    ;------------------------------------
 811    print(" 9. Plotting Psfc Regression coefficients")
 812    plts = new (2,"graphic")
 813
 814    res                      = True
 815    res@gsnMaximize          = True
 816    res@gsnDraw              = False
 817    res@gsnFrame             = False
 818    res@tmXTOn               = False
 819    res@tmYROn               = False
 820    res@cnFillOn             = Fill_Color
 821    res@lbLabelAutoStride     = True
 822    res@lbOrientation         = "Vertical"
 823    res@gsnSpreadColors      = True
 824    res@gsnSpreadColorStart  = 20
 825    res@gsnSpreadColorEnd    = -2
 826    res@gsnYAxisIrregular2Linear=True
 827
 828    if ( .not. By_Levels ) then
 829       res@tiYAxisString = "Sigma Values"
 830       res@trYReverse    = True
 831       res@tmYLMode      = "Explicit"
 832       res@tmYLValues    = sig_values_for_plot
 833       res@tmYLLabels    = sig_array_for_plot
 834    else
 835       res@tiYAxisString = "Sigma Levels"
 836       res@trYMinF       =  1.0
 837       res@trYMaxF       =  nsig + 0.0
 838    end if
 839    res@tiXAxisString = "Latitude"
 840    res@tmXBMode      = tmMode
 841    if ( res@tmXBMode .eq. "Explicit" ) then
 842       res@tmXBValues = lat_values_for_plot
 843       res@tmXBLabels = lat_array_for_plot
 844    end if
 845
 846    res@gsnLeftString = "10**9 x Psfc-Regression Coefficients"
 847    plts(0) = gsn_csm_contour (wks,reg_ps(:,1:nlat),res)
 848    delete (res)
 849
 850    ;--------------------------------------------------------------------
 851    ;  Plot Vertical profile of Psfc regression coeff
 852    ;--------------------------------------------------------------------
 853    reg_ps_avg = dim_avg_Wrap(reg_ps)
 854
 855    res                   = True
 856    res@gsnMaximize       = True
 857    res@gsnDraw           = False
 858    res@gsnFrame          = False
 859    res@tmXTOn            = False
 860    res@tmYROn            = False
 861    res@xyLineThicknesses = (/4.0,4.0/)
 862    res@xyDashPatterns    = (/0,4/)
 863    res@xyLineColors      = (/"blue"/)
 864
 865    res@tiXAxisString = "Regression Coefficient"
 866    if ( .not. By_Levels ) then
 867       res@tiYAxisString = "Sigma Values"
 868       res@trYReverse    = True
 869       res@tmYLMode      = "Explicit"
 870       res@tmYLValues    = sig_values_for_plot
 871       res@tmYLLabels    = sig_array_for_plot
 872    else
 873       res@tiYAxisString = "Sigma Levels"
 874       res@trYMinF       =  1.0
 875       res@trYMaxF       =  nsig + 0.0
 876    end if
 877
 878    res@xyExplicitLegendLabels = "ps"
 879    res@gsnLeftString   = "10**9 x Psfc - Averaged Regression Coefficients"
 880    if ( .not. By_Levels ) then
 881       plts(1) = gsn_csm_xy (wks,reg_ps_avg,rsig,res)
 882    else
 883       plts(1) = gsn_csm_xy (wks,reg_ps_avg,yc,res)
 884    end if
 885
 886    resP                     = True
 887    resP@gsnMaximize         = True
 888    resP@gsnPaperOrientation = "portrait"
 889    resP@gsnPanelYWhiteSpacePercent = 5.0  ; default 1.0
 890    resP@gsnPanelBottom      = 0.10
 891    gsn_panel(wks,plts,(/2,1/),resP)
 892
 893    delete (res)
 894    delete (resP)
 895    delete (plts)
 896
 897    ;------------------------------------------
 898    ; Plot Temperature Regression coefficients
 899    ;------------------------------------------
 900    print("10. Plotting Temperature Regression coefficients")
 901
 902    res                     = True
 903    res@cnFillOn            = Fill_Color
 904    res@gsnSpreadColors     = True
 905    res@gsnSpreadColorStart = 20
 906    res@gsnSpreadColorEnd   = -2
 907    res@lbLabelAutoStride   = True
 908    res@lbOrientation       = "Vertical"
 909    res@gsnDraw             = False
 910    res@gsnFrame            = False
 911    res@tmXTOn              = False
 912    res@tmYROn              = False
 913    res@gsnXAxisIrregular2Linear=True
 914    res@gsnYAxisIrregular2Linear=True
 915    if ( .not. By_Levels ) then
 916       ; Reverses X-axis
 917       res@tiXAxisString = "Sigma Values"
 918       res@tiYAxisString = "Sigma Values"
 919       res@trXReverse = True
 920       res@tmXBMode   = "Explicit"
 921       res@tmXBValues = sig_values_for_plot
 922       res@tmXBLabels = sig_array_for_plot
 923       ; Reverses Y-axis
 924       res@trYReverse = True
 925       res@tmYLMode   = "Explicit"
 926       res@tmYLValues = sig_values_for_plot
 927       res@tmYLLabels = sig_array_for_plot
 928    else
 929       res@tiXAxisString = "Sigma Levels"
 930       res@tiYAxisString = "Sigma Levels"
 931    end if
 932
 933    res@gsnLeftString = "10**9 x Averaged Temp Regression coeff "
 934    plts = gsn_csm_contour (wks,dim_avg_Wrap(reg_t),res) ; create plot
 935    draw(plts)
 936    frame(wks)
 937
 938    delete (res)
 939    delete(plts)
 940    destroy(wks)