<filename>components/elm/tools/clm4_5/OLMT/plotcase.py
#!/usr/bin/python

import os, sys, csv, glob
import numpy, scipy
from scipy.io import netcdf
from optparse import OptionParser
import matplotlib as mpl

def getvar(fname, varname, npf, index, scale_factor):
    usescipy = False
    try:
        import Scientific.IO.NetCDF as netcdf
    except ImportError:
        import scipy
        from scipy.io import netcdf
        usescipy = True
    if (usescipy):
        nffile = netcdf.netcdf_file(fname,"r",mmap=False)
        var = nffile.variables[varname]
        varvals = var[0:npf,index].copy() * scale_factor    #works for vector only?
        nffile.close()
    else:
        nffile = netcdf.NetCDFFile(fname,"r")
        var = nffile.variables[varname]
        varvals = var.getValue()[0:npf,index] * scale_factor
        nffile.close()
    return varvals


parser = OptionParser()
parser.add_option("--ad_Pinit", dest="ad_Pinit", default=False, action="store_true",\
                  help="Initialize AD spinup with P pools and use CNP mode")
parser.add_option("--csmdir", dest="mycsmdir", default='', \
                  help = 'Base CESM directory (default = ..)')
parser.add_option("--cases", dest="mycase", default='', \
                  help = "name of case id prefixs to plot (comma delmited)")
parser.add_option("--compset", dest="compset", default="I20TRCLM45CN", \
                  help = "Compset to plot")
parser.add_option("--titles", dest="titles", default='', \
                  help = "titles of case to plot (for legend)")
parser.add_option("--obs", action="store_true", default=False, \
                  help = "plot observations", dest="myobs")
parser.add_option("--sites", dest="site", default="none", \
                  help = 'site (to plot observations)')
parser.add_option("--timezone", dest="timezone", default=0, \
                  help = 'time zone (relative to UTC')
parser.add_option("--varfile", dest="myvarfile", default='varfile', \
                  help = 'file containing list of variables to plot')
parser.add_option("--vars", dest="myvar", default='', \
                  help="variable to plot (overrides varfile, " \
                  +"sends plot to screen")
parser.add_option("--avpd", dest="myavpd", default=1, \
                  help = 'averaging period in # of output timesteps' \
                  +' (default = 1)')
parser.add_option("--hist_mfilt", dest="myhist_mfilt", default=-999, \
                  help = 'beginning model year to plot')
parser.add_option("--hist_nhtfrq", dest="myhist_nhtfrq", default=-999, \
                  help = 'beginning model year to plot')
parser.add_option("--ystart", dest="myystart", default=1, \
                  help = 'beginning model year to plot')
parser.add_option("--yend", dest="myyend", default=9999, \
                  help = 'final model year to plot')
parser.add_option("--ystart_obs", dest="ystart_obs", default=0, \
                  help = 'beginning model year to plot')
parser.add_option("--yend_obs", dest="yend_obs", default=0, \
                  help = 'final model year to plot')
parser.add_option("--diurnal", dest="mydiurnal", default=False, \
                  action="store_true", help = 'plot diurnal cycle')
parser.add_option("--dstart", dest="dstart", default=1, \
                  help = 'beginning model DOY to plot (for diruanl average)')
parser.add_option("--dend", dest="dend", default=365, \
                  help = 'final model DOY to plot (for diurnal average)')
parser.add_option("--seasonal", dest="myseasonal", default=False, \
                  action="store_true", help = 'plot seasonal cycle')
parser.add_option("--h1", dest="h1", default=False, \
                  action="store_true", help = 'Use h1 history files')
parser.add_option("--h2", dest="h2", default=False, \
                  action="store_true", help = 'Use h2 history files')
parser.add_option("--h3", dest="h3", default=False, \
                  action="store_true", help = 'Use h3 history files')
parser.add_option("--h4", dest="h4", default=False, \
                  action="store_true", help = 'Use h4 history files')
parser.add_option("--index", dest="index", help = 'index (site or pft)', \
                   default=0)
parser.add_option("--spinup", dest="spinup", help = 'plot Ad and final spinup', \
                   default=False, action="store_true")
parser.add_option("--scale_factor", dest="scale_factor", help = 'scale factor', \
                   default=-999)
parser.add_option("--ylog", dest="ylog", help="log scale for Y axis", \
                  action="store_true", default=False)
parser.add_option("--pdf", dest="pdf", help="save plot to pdf", \
                  action="store_true", default=False)
parser.add_option("--noplot", dest="noplot", help="Do not make plots", \
                  action="store_true", default=False)
(options,args) = parser.parse_args()
               

cesmdir=os.path.abspath(options.mycsmdir)                 

if (options.pdf):
    mpl.use('Agg')	
import matplotlib.pyplot as plt

mycases = options.mycase.split(',')
mysites = options.site.split(',')
mycompsets = options.compset.split(',')

ncases = 1
if (len(mycases) > 1):
  ncases = len(mycases)
  mysites=[]
  for c in range(0,ncases):
    mysites.append(options.site)
    #if (len(mycompsets) == 1):
    mycompsets.append(options.compset)
  mytitles = mycases
elif (len(mysites) > 1):
  ncases = len(mysites)
  mycases=[]
  mycompsets=[]
  for c in range(0,ncases):
    mycases.append(options.mycase)
    mycompsets.append(options.compset)
  mytitles = mysites
elif (len(mycompsets) > 1):
  ncases = len(mycompsets)
  mycases=[]
  mysites=[]
  for c in range(0,ncases):
    mycases.append(options.mycase)
    mysites.append(options.site)  
  mytitles = mycompsets
else:
  mytitles=[]
  mytitles.append(mysites[0])

if (options.titles != ''):
  mytitles = options.titles.split(',')

obs     = options.myobs
myobsdir = '/lustre/or-hydra/cades-ccsi/scratch/dmricciuto/fluxnet'

#get list of variables from varfile
myvars=[]

if (options.myvar == ''):
    if os.path.isfile('./'+options.myvarfile):
        input = open('./'+options.myvarfile)
        for s in input:
            myvars.append(s.strip())
    else:
        print('Error:  invalid varfile')
        sys.exit()
    terminal = 'postscript'
else:
    terminal=''
    myvars = options.myvar.split(',')
    

avpd      = int(options.myavpd)        # desired averaging period in output timestep
ystart    = int(options.myystart)      # beginning year to plot/average
yend      = int(options.myyend)        # final year to plot/average 
yend_all  = yend                       # keep track of last year for which datasets exist
mylat_vals =[]
mylon_vals= []

avtype = 'default'
if (options.mydiurnal):
    avtype = 'diurnal'
    avpd=1
if (options.myseasonal):
    avtype = 'seasonal'

#------------------------------------------------------------------------------

#site = options.site
#compset = options.compset

#dirs=[]
nvar = len(myvars)    
x_toplot    = numpy.zeros([ncases, 2000000], numpy.float)
data_toplot = numpy.zeros([ncases, nvar, 2000000], numpy.float)
obs_toplot  = numpy.zeros([ncases, nvar, 2000000], numpy.float)+numpy.NaN
err_toplot  = numpy.zeros([ncases, nvar, 2000000], numpy.float)+numpy.NaN
snum        = numpy.zeros([ncases], numpy.int)

for c in range(0,ncases):
    if (mycases[c] == ''):
        mydir = cesmdir+'/'+mysites[c]+'_'+mycompsets[c]+'/run/'
    else:
        mydir = cesmdir+'/'+mycases[c]+'_'+mysites[c]+'_'+mycompsets[c]+'/run/'
    print('Processing '+mydir)

    #query lnd_in file for output file information
    if ((options.myhist_mfilt == -999 or options.myhist_nhtfrq == -999)):
        #print('Obtaining output resolution information from lnd_in')
        input = open(mydir+"/lnd_in")
        npf=-999
        tstep=-999
        input = open(mydir+"/lnd_in")
        for s in input:
	    if ('hist_mfilt' in s):
	        mfiltinfo = s.split()[2]
                npf = int(mfiltinfo.split(',')[0])
                if (options.h1): 
                    npf = int(mfiltinfo.split(',')[1])
                if (options.h2):
                    npf = int(mfiltinfo.split(',')[2])
                if (options.h3):
                    npf = int(mfiltinfo.split(',')[3])
                if (options.h4):
                    npf = int(mfiltinfo.split(',')[4])
            if ('hist_nhtfrq' in s):
                nhtfrqinfo = s.split()[2]
                tstep = int(nhtfrqinfo.split(',')[0])
                if (options.h1):
                    tstep = int(nhtfrqinfo.split(',')[1])
                if (options.h2):
                    tstep = int(nhtfrqinfo.split(',')[2])
                if (options.h3):
                    tstep = int(nhtfrqinfo.split(',')[3])
                if (options.h4):
                    tstep = int(nhtfrqinfo.split(',')[4])
        input.close()
    else:
        npf   = int(options.myhist_mfilt)
        tstep = int(options.myhist_nhtfrq)
   
    if (npf == -999 or tstep == -999):
        print('Unable to obtain output file information from lnd_in.  Exiting')
        sys.exit()

    yststr=str(100000+ystart)
    #determine type of file to plot
    if (options.h4):
        hst = 'h4'
    elif (options.h3):
        hst = 'h3'
    elif (options.h2):
        hst = 'h2'
    elif (options.h1):
        hst = 'h1'
    else:
        hst = 'h0'

    if (tstep == 0):
        ftype = 'default'
        mytstep = 'monthly'
        npy=12
    else:
        ftype = 'custom'
        if (abs(npf) == 8760):
            mytstep = 'halfhourly'
            npy=8760
        elif (abs(npf) == 365):
            mytstep = 'daily'
            npy=365
        elif (abs(npf) == 1):
            mytstep = 'annual'
            npy=1
        nhtot=-1*tstep*npf
        nypf = max(1, nhtot/8760)

    #initialize data arrays
    mydata      = numpy.zeros([nvar,2000000], numpy.float)
    myobs       = numpy.zeros([nvar,2000000], numpy.float)+numpy.NaN
    myerr       = numpy.zeros([nvar,2000000], numpy.float)+numpy.NaN
    x           = numpy.zeros([2000000], numpy.float)
    nsteps=0
 
    if (c == 0):   
        var_units=[]
        var_long_names=[]
        myscalefactors=[]

    #Get observations
    if (obs):
        myobsfiles = os.listdir(myobsdir+'/'+mytstep+'/')
        for f in myobsfiles:
            if mysites[c] in f and '.csv' in f:
                myobsfile = myobsdir+'/'+mytstep+'/'+f
        avpd_obs = 1
        if (mytstep == 'halfhourly' and '_HH_' in myobsfile):
            avpd_obs = 2 
        if (os.path.exists(myobsfile) and ystart < 1900):
            print ('Getting start and end year information from observation file')
            thisrow=0
            myobs_input = open(myobsfile)
            for s in myobs_input:
                if thisrow == 1:
                    ystart = int(s[0:4])
                elif (thisrow > 1):
                    if (int(options.myyend) > 9000):
                        yend = min(int(s[0:4]), int(options.myyend))
                thisrow=thisrow+1
            myobs_input.close
        for v in range(0,nvar):
            if (os.path.exists(myobsfile)):
                myobs_in = open(myobsfile)
                thisrow=0
                thisob=0
                for s in myobs_in:
                    if (thisrow == 0):
                        header = s.split(',')
                    else:
                        myvals = s.split(',')
                        thiscol=0
                        if int(myvals[0][0:4]) >= ystart and int(myvals[0][0:4]) <= yend:
                            if (thisob == 0):
                                thisob = (int(myvals[0][0:4])-ystart)*npy*avpd_obs
                            if (thisob % avpd_obs ==  0):
                                myobs[v,thisob/avpd_obs] = 0.0
                                myerr[v,thisob/avpd_obs] = 0.0
                            for h in header:
                                if (h.strip() == 'NEE_CUT_REF' and 'NEE' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] + float(myvals[thiscol])/avpd_obs
                                elif (h.strip () == 'NEE_CUT_REF_JOINTUNC' and \
                                      'NEE' in myvars[v]):
                                    myerr[v,thisob/avpd_obs] = myerr[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'GPP_NT_CUT_REF' and 'GPP' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'GPP_NT_CUT_SE' and 'GPP' in myvars[v]):
                                    myerr[v,thisob/avpd_obs] = myerr[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'RECO_NT_CUT_REF' and 'ER' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'RECO_NET_CUT_SE' and 'ER' in myvars[v]):
                                    myerr[v,thisob/avpd_obs] = myerr[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'LE_F_MDS' and 'EFLX_LH_TOT' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'LE_RANDUNC' and 'EFLX_LH_TOT' in myvars[v]):
                                    myerr[v,thisob/avpd_obs] = myerr[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'H_F_MDS' and 'FSH' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'H_RANDUNC' and 'FSH' in myvars[v]):
                                    myerr[v,thisob/avpd_obs] = myerr[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'TA_F_MDS' and 'TBOT' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'SWIN_F_MDS' and 'FSDS' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'WS_F' and 'WIND' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                elif (h.strip() == 'P_F' and 'RAIN' in myvars[v]):
                                    myobs[v,thisob/avpd_obs] = myobs[v,thisob/avpd_obs] +float(myvals[thiscol])/avpd_obs
                                if myobs[v,thisob/avpd_obs] < -4000:
                                    myobs[v,thisob/avpd_obs] = numpy.NaN
                                thiscol=thiscol+1
                            thisob=thisob+1
                    thisrow = thisrow+1
                myobs_in.close()
            else:
                print 'Error reading observations for '+mysites[c]

    #read monthly .nc files (default output)
    if (ftype == 'default'):
        for v in range(0,nvar):
            nsteps = 0
            for y in range(ystart,yend+1):
                yst=str(10000+y)[1:5]
                for m in range(0,12):
                    mst=str(101+m)[1:3]
                    myfile = os.path.abspath(mydir+'/'+mycases[c]+'_'+mysites[c]+'_'+mycompsets[c]+ \
                                             ".clm2."+hst+"."+yst+"-"+mst+".nc")
                    #get units/long names from first file
                    if (os.path.exists(myfile)):
                        if (y == ystart and m == 0 and c == 0):
                            nffile = netcdf.netcdf_file(myfile,"r")
                            varout=nffile.variables[myvars[v]]
                            var_long_names.append(varout.long_name)
                            nffile.close()
                            if (float(options.scale_factor) < -900):
                                if ('gC/m^2/s' in varout.units):
                                    myscalefactors.append(3600*24)
                                    var_units.append('g.C/m2/day')
                                else:
                                    myscalefactors.append(1.0)
                                    var_units.append(varout.units.replace('^',''))
                            else:
                                myscalefactors.append(float(options.scale_factor))
                                var_units.append(varout.units.replace('^',''))
                        
                        if (y == ystart and m == 0 and v == 0):      # get lat/lon info
                            nffile = netcdf.netcdf_file(myfile,"r")
                            mylat_vals.append(nffile.variables['lat'][0])
                            mylon_vals.append(nffile.variables['lon'][0])
                            nffile.close()
                        x[nsteps] = y+m/12.0
                        myvar_temp = getvar(myfile, myvars[v],npf,int(options.index), \
                                            myscalefactors[v])
                        mydata[v,nsteps] = myvar_temp
                        if (myvars[v] == 'RAIN'):
                            myvar_temp2 = getvar(myfile,'SNOW',npf,int(options.index), \
                                            myscalefactors[v])
                            mydata[v,nsteps] = mydata[v,nsteps]+myvar_temp2
                    else:
                        if (v == 0 and m == 0):
                           print 'Warning: '+myfile+' does not exist'
                        x[nsteps] = y+m/12.0
                        mydata[v,nsteps] = numpy.NaN
                        if (y-1 < yend_all):
                            yend_all = y-1
                    nsteps = nsteps+1

    #read annual .nc files
    if (ftype == 'custom'):
        for v in range(0,nvar):
            nsteps=0
            nfiles = (yend-ystart)/nypf
            nc=1
            starti=0
            ylast=0
            if (options.spinup):
                nc=2
            if (npf == 1):
                starti = 1
            for n in range(0,nc):
                if ((options.spinup)and n== 0):
                    if (mycases[c] == ''):
                        if (options.ad_Pinit):
                            mydir = cesmdir+'/'+mysites[c]+'_'+mycompsets[c]+'_ad_spinup/run/'
                        else:
                            mydir = cesmdir+'/'+mysites[c]+'_'+mycompsets[c].replace('CNP','CN')+ \
	                          '_ad_spinup/run/'
                    else:
                        if (options.ad_Pinit):
                            mydir = cesmdir+'/'+mycases[c]+'_'+mysites[c]+'_'+ \
                                    mycompsets[c]+'_ad_spinup/run/'
                            thiscompset = mycompsets[c]+'_ad_spinup'
                        else:
                            mydir = cesmdir+'/'+mycases[c]+'_'+mysites[c]+'_'+ \
                                    mycompsets[c].replace('CNP','CN')+'_ad_spinup/run/'
                            thiscompset = mycompsets[c].replace('CNP','CN')+'_ad_spinup'
                else:
                    if (mycases[c] == ''):
                        mydir = cesmdir+'/'+mysites[c]+'_'+mycompsets[c]+'/run/'
                    else:
                        mydir = cesmdir+'/'+mycases[c]+'_'+mysites[c]+'_'+ \
                                mycompsets[c]+'/run/'
                    thiscompset = mycompsets[c]
                for y in range(starti,nfiles+1):   #skip first file in spinup
                    yst=str(10000+ystart+(y*nypf))[1:5]
                    if (mycases[c].strip() == ''):
                        myfile = os.path.abspath(mydir+'/'+mycases[c]+'_'+thiscompset+".clm2."+hst+ \
                                                 "."+yst+"-01-01-00000.nc")
                    else:
                        myfile = os.path.abspath(mydir+'/'+mycases[c]+"_"+mysites[c]+'_'+thiscompset+ \
                                                 ".clm2."+hst+"."+yst+"-01-01-00000.nc")
                    if (os.path.exists(myfile)):
                        if (n == 0):
                            ylast = y
                        if (y == starti and n == 0 and c == 0):
                            nffile = netcdf.netcdf_file(myfile,"r")
                            varout=nffile.variables[myvars[v]]
                            var_long_names.append(varout.long_name)
                            if (float(options.scale_factor) < -900):
                                if ('gC/m^2/s' in varout.units):
                                    if (npf >= 365):
                                        myscalefactors.append(3600*24)
                                        var_units.append('g.C/m2/day')
                                    else:
                                        myscalefactors.append(3600*24*365)
                                        var_units.append('g.C/m2/yr')
                                else:
                                    myscalefactors.append(1.0)
                                    var_units.append(varout.units.replace('^',''))
                            else:
                                 myscalefactors.append(float(options.scale_factor))
                                 var_units.append(varout.units.replace('^',''))
                            nffile.close()
                        if (y == starti and n == 0 and v == 0):      # get lat/lon info
                            nffile = netcdf.netcdf_file(myfile,"r")
                            mylat_vals.append(nffile.variables['lat'][0])
                            mylon_vals.append(nffile.variables['lon'][0])
                            nffile.close()

                        myvar_temp = getvar(myfile,myvars[v],npf,int(options.index), \
                                            myscalefactors[v])
                        if (myvars[v] == 'RAIN'):
                            myvar_temp2 = getvar(myfile,'SNOW',npf,int(options.index), \
                                            myscalefactors[v])

                        if len(myvar_temp) == npf:
                            for i in range(0,npf):
                                myind = ylast*n*npf+y*npf+i
                                x[nsteps] = ystart+(ylast*n*nypf+y*nypf) + nypf*(i*1.0-0.5)/npf
                                mydata[v,nsteps] = myvar_temp[i]
                                if (myvars[v] == 'RAIN'):    #add snow for total precip
                                    mydata[v,nsteps] = mydata[v,nsteps]+myvar_temp2[i]
                                nsteps=nsteps+1
                        else:
                            for i in range(0,npf):
                                myind=ylast*n*npf+(y-1)*npf+i
                                x[myind] = ystart+(ylast*n*nypf+y*nypf) + nypf*(i*1.0-0.5)/npf
                                mydata[v,myind] = numpy.NaN
                                nsteps=nsteps+1 
  		    else:	
                         if (v == 0):
                             print('Warning: '+myfile+' does not exist')
                         if (y-1 < yend_all):
                             yend_all = y-1
                         for i in range(0,npf):
                             if (n == nc-1):
                                 myind=ylast*n*npf+y*npf+i
                                 x[myind] = ystart+(ylast*n*nypf+y*nypf) + nypf*(i*1.0-0.5)/npf
                                 mydata[v,myind] = numpy.NaN
                                 nsteps=nsteps+1

    #perform averaging and write output files 
    if (avtype == 'default'):
        for v in range(0,nvar):
            snum[c] = 0
            for s in range(0,int(nsteps/avpd)): 
                x_toplot[c, snum[c]]       = sum(x[s*avpd:(s+1)*avpd])/avpd
                data_toplot[c, v, snum[c]] = sum(mydata[v,s*avpd:(s+1)*avpd])/avpd
                obs_toplot[c, v, snum[c]]  = sum(myobs[v,s*avpd:(s+1)*avpd])/avpd
                if (min(myerr[v,s*avpd:(s+1)*avpd]) < -9000):
                    err_toplot[c,v,snum[c]] = 0
                else:
                    err_toplot[c, v, snum[c]]  = sum(myerr[v,s*avpd:(s+1)*avpd])/avpd
                snum[c] = snum[c]+1
          

    #diurnal average (must have hourly output)
    if (avtype == 'diurnal'):
        snum[c]=24
        for v in range(0,nvar):
            mysum = numpy.zeros(snum[c], numpy.float)
            mysum_obs = numpy.zeros(snum[c], numpy.float)
            myct = numpy.zeros(snum[c],numpy.float)
            myct_obs = numpy.zeros(snum[c],numpy.float)
            for y in range(0,(yend_all-ystart+1)):
                for d in range (int(options.dstart),int(options.dend)):
                    for s in range(0,snum[c]):        
                        h=s
                        if (h >= 24):
                            h=h-24
                        mysum[s] = mysum[s]+mydata[v,y*8760+(d-1)*24+h-int(options.timezone)+1]
                        myct[s] = myct[s]+1
                        if (myobs[v,y*8760+(d-1)*24+h] > -900):
                            mysum_obs[s] = mysum_obs[s]+myobs[v,y*8760+(d-1)*24+h]
                            myct_obs[s] = myct_obs[s]+1
            for s in range(0,snum[c]):
                if (myct_obs[s] > 0):
                    mysum_obs[s] = mysum_obs[s]/myct_obs[s]
                else:
                    mysum_obs[s] = numpy.NaN
                x_toplot[c,s] = s+1
                obs_toplot[c, v, s] = mysum_obs[s]
                data_toplot[c, v, s] = mysum[s]/myct[s]

    #seasonal average (assumes default monthly output)
    if (avtype == 'seasonal'):
        for v in range(0,nvar):
            snum[c] = 12
            mysum=numpy.zeros(snum[c], numpy.float)
            mysum_obs = numpy.zeros(snum[c], numpy.float)
            mycount_obs = numpy.zeros(snum[c], numpy.int)
            for y in range(0,(yend_all-ystart+1)):
                for s in range(0,snum[c]):
                    mysum[s]=mysum[s]+mydata[v,(y*12+s)]/float(yend_all-ystart+1)
                    if (myobs[v,(y*12+s)] > -900):
                        mysum_obs[s]=mysum_obs[s]+myobs[v,(y*12+s)]
                        mycount_obs[s] = mycount_obs[s]+1
            for s in range(0,snum[c]):
                if (mycount_obs[s] > 0):
                    mysum_obs[s] = mysum_obs[s]/mycount_obs[s]
                else:
                    mysum_obs[s] = numpy.NaN
                x_toplot[c,s] = s+1.5
                obs_toplot[c,v,s]  = mysum_obs[s]
                data_toplot[c,v,s] = mysum[s]

                
#diagnostics, outputs and plots

if (options.spinup):
    analysis_type = 'spinup'
elif (options.mydiurnal):
    analysis_type = 'diurnalcycle_'+str(options.dstart)+'_'+str(options.dend)
elif (options.myseasonal):
    analysis_type = 'seasonalcycle'
elif (mytstep == 'halfhourly'):
    analysis_type = 'hourly'
else:
    analysis_type=mytstep

rmse = numpy.zeros([len(myvars),ncases],numpy.float)
bias = numpy.zeros([len(myvars),ncases],numpy.float)
corr = numpy.zeros([len(myvars),ncases],numpy.float)

for v in range(0,len(myvars)):
    if (not options.noplot):
        fig = plt.figure()
        ax = plt.subplot(111)
    colors=['b','g','r','c','m','y','k','b','g','r','c','m','y','k','b','g','r','c','m','y','k']
    styles=['-','-','-','-','-','-','-','--','--','--','--','--','--','--','-.','-.','-.','-.','-.','-.','-.']
    for c in range(0,ncases):
        #Output data in netcdf format
        if (c == 0):
            if (v == 0):
                ftype_suffix=['model','obs']
                os.system('mkdir -p ./plots/'+mycases[0]+'/'+analysis_type)
                for ftype in range(0,2):
                    outdata = netcdf.netcdf_file('./plots/'+mycases[0]+'/'+analysis_type+'/'+mycases[0]+"_"+mysites[0]+'_'+ftype_suffix[ftype]+".nc","w",mmap=False)
                    outdata.createDimension('time',snum[c])
                    #outdata.createDimension('lat',ncases)
                    #outdata.createDimension('lon',ncases)
                    outdata.createDimension('gridcell',ncases)
                    outdata.createDimension('strlen',6)
                    mylat = outdata.createVariable('lat','f',('gridcell',))
                    mylat.long_name='coordinate latitude'
                    mylat.units='degrees_north'
                    mylon = outdata.createVariable('lon','f',('gridcell',))
                    mylon.long_name='coordinate longitude'
                    mylon.units='degrees_east'
                    mytime = outdata.createVariable('time','f',('time',))
                    mytime.long_name='time'
                    mytime.units='days since '+str(ystart)+'-01-01 00:00:00'
                    mytime.calendar='noleap'
                    mytime[:] = (x_toplot[0,0:snum[c]]-ystart)*365
                    #myname = outdata.createVariable('site_name','c',('lat','lon','strlen'))
                    myname = outdata.createVariable('site_name','c',('gridcell','strlen'))
                    myname[:,:] = ''   #changed for gridcell
                    outdata.close()
        for ftype in range(0,2):
            outdata = netcdf.netcdf_file('./plots/'+mycases[0]+'/'+analysis_type+'/'+mycases[0]+"_"+mysites[0]+'_'+ftype_suffix[ftype]+".nc","a",mmap=False)
            if (c == 0):
                #myvar = outdata.createVariable(myvars[v],'f',('time','lat','lon'))
                myvar = outdata.createVariable(myvars[v],'f',('time','gridcell'))
                myvar.units=var_units[v]
                myvar.missing_value=1e36
                myvar[:,:]=myvar.missing_value   #changed for gridcell
            else:
                myvar=outdata.variables[myvars[v]]
            scalefac = 1.0
            if (var_units[v] == 'g.C/m2/day'):
                myvar.units = 'kg.C/m2/s'
                scalefac = 1.0 / (3600*24*1000.0)
            if (ftype == 0):
                myvar[:,c] = data_toplot[c,v,0:snum[c]]*scalefac  #changed for gridcell
            if (ftype == 1):
                myvar[:,c] = obs_toplot[c,v,0:snum[c]]*scalefac   #changed for gridcell
            if (v == 0):
                myname = outdata.variables['site_name']
                myname[c,0:6] = str(mysites[c])[0:6]   #changed for gridcell
                mylat = outdata.variables['lat']
                mylat[c] = mylat_vals[c]
                mylon = outdata.variables['lon']
                mylon[c] = mylon_vals[c]
            outdata.close()
        
        #----------------------- plotting ---------------------------------
        if (options.noplot == False):
            gind=[]
            for i in range(0,snum[c]):
                if (obs_toplot[c,v,i] < -900):
                    obs_toplot[c,v,i] = numpy.nan
                else:
                    gind.append(i)
                    rmse[v,c] = rmse[v,c] + (data_toplot[c,v,i]-obs_toplot[c,v,i])**2.0
                    bias[v,c] = bias[v,c] + (data_toplot[c,v,i]-obs_toplot[c,v,i])
            rmse[v,c] = (rmse[v,c]/len(gind))**0.5
            bias[v,c] = bias[v,c]/len(gind)
            corr[v,c] = numpy.corrcoef(data_toplot[c,v,gind],obs_toplot[c,v,gind])[1,0]

            if (options.ylog):
                ax.plot(x_toplot[c, 0:snum[c]], abs(data_toplot[c,v,0:snum[c]]), label=mytitles[c], color=colors[c], \
                  linestyle=styles[c], linewidth=3)
            else:
                ax.plot(x_toplot[c, 0:snum[c]], (data_toplot[c,v,0:snum[c]]), label=mytitles[c], color=colors[c], \
	          linestyle=styles[c], linewidth=3)
                if (options.myobs and c == 0):
                    ax.errorbar(x_toplot[c, 0:snum[c]], obs_toplot[c,v,0:snum[c]], yerr = err_toplot[c,v,0:snum[c]], \
                                color=colors[c], fmt='o')

    if (options.noplot == False):
        if (avtype == 'seasonal'):
            plt.xlabel('Model Month')
        elif (avtype == 'diurnal'):
            plt.xlabel('Model Hour (LST)')
        else:
            plt.xlabel('Model Year')
 
        plt.ylabel(myvars[v]+' ('+var_units[v]+')')
        box = ax.get_position()
        ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5),prop={'size': 10})
        plt.title(var_long_names[v]+' at '+mysites[0])
        if (options.ylog):
            plt.yscale('log')

        os.system('mkdir -p ./plots/'+mycases[0]+'/'+analysis_type)
        if (options.pdf):
            fig_filename = './plots/'+mycases[0]+'/'+analysis_type+'/'+mysites[0]+'_'+myvars[v]
            fig_filename = fig_filename+'_'+analysis_type
            fig.savefig(fig_filename+'.pdf')

if (not options.pdf and not options.noplot):
    plt.show()
