source: trunk/UTIL/PYTHON/gcm.py @ 389

#!/usr/bin/env python

### T. Navarro + A. Spiga + A. Colaitis

###########################################################################################
###########################################################################################
### What is below relate to running the file as a command line executable (very convenient)
if __name__ == "__main__":
    import sys
    from optparse import OptionParser    ### to be replaced by argparse
    from api_wrapper import api_onelevel
    from zrecast_wrapper import call_zrecast
    from netCDF4 import Dataset
    from myplot import getlschar, separatenames, readslices, adjust_length
    from os import system
    from planetoplot import planetoplot
    from myscript import getparseroptions
    import numpy as np

    #############################
    ### Get options and variables
    parser = OptionParser()
    getparseroptions(parser)
    (opt,args) = parser.parse_args()

    if opt.file is None: 
        print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
        exit()
    if opt.var is None and opt.anomaly is True:
        print "Cannot ask to compute anomaly if no variable is set"
        exit()
    print "Options:", opt

#    listvar = ''
#    if opt.var is None:
#        zerange = [-999999]
#    else:
#        zelen = len(opt.var)
#        zerange = range(zelen)
#        #if zelen == 1: listvar = opt.var[0] + ','
#        #else         :
#        for jjj in zerange: listvar += opt.var[jjj] + ','
#        listvar = listvar[0:len(listvar)-1]
#        vmintab = adjust_length (opt.vmin, zelen)
#        vmaxtab = adjust_length (opt.vmax, zelen)


    ################################
    ### General check
 
    if opt.fref is not None:
       if opt.operat is None:
          print "you must specify an operation when using a reference file"
          exit()
    if opt.operat in ["+","-"]:
       if opt.fref is None:
          print "you must specifiy a reference file when using inter-file operations"
          exit()

    interpref=False
    if opt.fref is not None:
       if opt.operat is not None:
          if opt.itp is not None:
             interpref=True

    ################################


    print "file, length", opt.file, len(opt.file)

    zeslat  = readslices(opt.slat)
    zeslon  = readslices(opt.slon)
    zesvert = readslices(opt.svert)
    zestime = readslices(opt.stime)
    print "slat,zeslat", opt.slat, zeslat
    print "slon,zeslon", opt.slon, zeslon
    print "svert,zesvert", opt.svert, zesvert
    print "stime,zestime", opt.stime, zestime

    for i in range(len(opt.file)):
      for j in range(len(opt.var)):

        zefiles = separatenames(opt.file[i])
        print "zefiles", zefiles
        
        if opt.vmin is not None : zevmin  = opt.vmin[min(i,len(opt.vmin)-1)]
        else: zevmin = None
        if opt.vmax is not None : zevmax  = opt.vmax[min(i,len(opt.vmax)-1)]
        else: zevmax = None
        print "vmin, zevmin", opt.vmin, zevmin
        print "vmax, zevmax", opt.vmax, zevmax
        
        zevar = separatenames(opt.var[j])
        zevars = zevar[j]
        zevar = zevar[0]
        print "var, zevar", opt.var, zevar
        
        #checkcoherence(len(zefiles),len(opt.slat),len(opt.slon),len(opt.stime))
        
        zefile = zefiles[0]
              
        #zelevel = opt.lvl   
        stralt = None
        #[lschar,zehour,zehourin] = getlschar ( zefile )  ## getlschar from wrfout (or simply return "" if another file)
        [lschar,zehour,zehourin] = ["",0,0]  ## dummy
    
        inputnvert = separatenames(opt.lvl)
        if np.array(inputnvert).size == 1:
            zelevel = float(inputnvert[0])
            ze_interp_levels = [-9999.]
        else:
            zelevel = -99.
            ze_interp_levels = np.linspace(float(inputnvert[0]),float(inputnvert[1]),float(inputnvert[2]))
        print 'level: ', zelevel
        print 'interp_levels: ',ze_interp_levels

        #####################################################
        ### Call Fortran routines for vertical interpolations  -->  zrecast for GCM ?       
#        if opt.itp is not None:
#            if zelevel == 0. and opt.itp == 4:  zelevel = 0.010
#            ### winds or no winds
#            if opt.winds            :  zefields = 'uvmet'
#            else                    :  zefields = ''
#            ### var or no var
#            #if opt.var is None      :  pass
#            if zefields == ''       :  zefields = listvar 
#            else                    :  zefields = zefields + "," + listvar 
#            if opt.var2 is not None : zefields = zefields + "," + opt.var2  
#            print zefields
#            zefile = api_onelevel (  path_to_input   = '', \
#                                     input_name      = zefile, \
#                                     fields          = zefields, \
#                                     interp_method   = opt.itp, \
#                                     interp_level    = ze_interp_levels, \
#                                     onelevel        = zelevel, \
#                                     nocall          = opt.nocall )
#            print zefile
#            zelevel = 0 ## so that zelevel could play again the role of nvert

# A.C.  #####################################################
        ### Call Fortran routines for vertical interpolations  -->  zrecast

        if opt.itp is not None:
           interpolated_files=""
           interpolated_files=call_zrecast(interp_mode=opt.itp,\
                    input_name=zefiles,\
                    fields=zevars)

           zefiles=interpolated_files
        if interpref:
           interpolated_ref=""
           interpolated_ref=call_zrecast(interp_mode=opt.itp,\
                    input_name=[opt.fref],\
                    fields=zevars)

           reffile=interpolated_ref[0]
        else:
           reffile=opt.fref
# Divers ####################################################
    
        zexaxis=[opt.xmin,opt.xmax]
        zeyaxis=[opt.ymin,opt.ymax]

        #############
        ### Main call
        name = planetoplot (zefiles,level=int(zelevel),vertmode=opt.itp,\
                proj=opt.proj,back=opt.back,target=opt.tgt,stride=opt.ste,var=zevar,\
                numplot=opt.num,colorb=opt.clb,winds=opt.winds,\
                addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
                itstep=opt.it,hole=opt.hole,save=opt.save,\
                anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.frt,\
                mult=opt.mult,zetitle=opt.zetitle,\
                slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime,\
                outputname=opt.out,resolution=opt.res,\
                ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
                invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
                blat=opt.blat,tsat=opt.tsat)
        print 'Done: '+name
        system("rm -f to_be_erased")
  
    #########################################################
    ### Generate a .sh file with the used command saved in it
    command = ""
    for arg in sys.argv: command = command + arg + ' '
    name = 'pycommand'
    f = open(name+'.sh', 'w')
    f.write(command)