Changeset 310 for trunk/MESOSCALE


Ignore:
Timestamp:
Sep 29, 2011, 3:59:01 PM (13 years ago)
Author:
aslmd
Message:

MESOSCALE:
-- Code in storm scenario corresponds to 'OMEGA' reference case [Julien Faure]
-- Easier settings for dust lifting without the need to recompile [see below]
-- A few modifications to plot and dust-devil detection PYTHON routines

29/09/11 == AS

--> To easily explore sensitivity to lifting thresholds: in dustlift.F, ustar_seuil=sqrt(stress/rho)

and alpha_lift[dust_mass] can be prescribed through an external stress.def parameter file.
--- alpha_lift[dust_number] is computed from alpha_lift[dust_mass] as in initracer.F
--- ustar_seuil is more user-friendly than stress because direct comparison with ustar from model

--> For the moment this is MESOSCALE only, but potentially useful to everyone

Location:
trunk/MESOSCALE/LMD_MM_MARS/SRC
Files:
2 added
1 deleted
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/2dplot.py

    r296 r310  
    2626#title = "Atmospheric temperature"
    2727
     28name = "wrfout_d01_2024-03-22_01:00:00_z"
     29itime = 12
     30ndiv = 10
     31zey = 120
     32var = "VMR_ICE"
     33title = "Volume mixing ratio of water ice [ppm]"
     34vmin = -0.5
     35vmax = 400.
     36
    2837nc = Dataset(name)
    2938
    3039what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d2=zey )
     40
    3141
    3242y = nc.variables["vert"][:]
    … …  
    3444horinp = len(what_I_plot[0,:])
    3545x = np.linspace(0.,horinp*500.,horinp) / 1000.
     46xlabel("Horizontal coordinate (km)")
     47
     48horinp = len(what_I_plot[0,:])
     49x = np.linspace(0.,horinp,horinp)
     50xlabel("x grid points")
     51
    3652
    3753zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
    … …  
    4662                       extend='neither',spacing='proportional')
    4763ptitle(title)
    48 xlabel("Horizontal coordinate (km)")
    4964ylabel("Altitude (m)")
    5065makeplotres(var+str(itime),res=200.,disp=False)

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py

    r251 r310  
    2323    from scipy.ndimage.measurements import minimum_position
    2424    from netCDF4 import Dataset
    25    
     25    import matplotlib.pyplot as plt
     26    import myplot as myp
     27   
    2628    ### LOAD NETCDF DATA
    2729    filename = "psfc.nc"
    … …  
    3436    allsizesx = []
    3537    allsizesy = []
    36     for i in range(0,shape[0]):
    37     #for i in range(0,2):
    38    
    39         print i, ' on ', shape[0]
     38    stride = 1 #5
     39    for i in range(0,shape[0],stride):
     40
    4041        psfc2d = np.array ( psfc [ i, : , : ] )
    4142   
    4243        ############### CRITERION
    43         ave = np.mean(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy
    44         where = np.where(psfc2d - ave < -0.3)  ## comme le papier Phoenix
     44        ave = np.mean(psfc2d,dtype=np.float64)  ## dtype otherwise inaccuracy
     45
     46        #limdp = -0.2 ## on en loupe pas mal
     47        #where = np.where(psfc2d - ave < limdp)  ## comme le papier Phoenix
    4548   
    46         std = np.std(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy
    47         fac = 2.5  ## not too low, otherwise vortices are not caught. 4 good choice.
     49        std = np.std(psfc2d,dtype=np.float64)   ## dtype otherwise inaccuracy
     50        fac = 4.   ## how many sigmas. not too low, otherwise vortices are not caught. 4 good choice.
     51                   ## 2.5 clearly too low, 3.5 not too bad, 4 probably good
     52        fac = 3.5
     53        #fac = 2.5
    4854        lim = ave - fac*std
    49         print lim, ave, std
    5055        where = np.where(psfc2d < lim)
    5156        ############### END CRITERION
    52     
     57   
    5358        lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine
    5459        lab[where] = 1.  ## do not treat points close to 'mean' (background) pressure
    55    
     60
     61        draw = False
     62        if draw:
     63        ##################################################################################
     64            vmin = -0.3
     65            vmax =  0.0
     66            ndiv = 3
     67            palette = plt.get_cmap(name="YlGnBu")
     68            what_I_plot = psfc2d-ave
     69            zevmin, zevmax = myp.calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
     70            what_I_plot = myp.bounds(what_I_plot,zevmin,zevmax)
     71            zelevels = np.linspace(zevmin,zevmax)
     72            fig = plt.figure()
     73            sub = myp.definesubplot(2,fig)
     74            plt.subplot(sub)
     75            plt.contourf(what_I_plot,zelevels,cmap=palette)
     76            plt.colorbar(fraction=0.05,pad=0.03,format="%.1f",\
     77                           ticks=np.linspace(zevmin,zevmax,ndiv+1),\
     78                           extend='both',spacing='proportional')
     79            plt.subplot(sub+1)
     80            palette = plt.get_cmap(name="binary")
     81            plt.contourf(lab,2,cmap=palette)
     82            plt.show()
     83        ##################################################################################
     84 
    5685        xx = []
    5786        yy = []
    … …  
    6998        sizex = detsize( xx, res = 10, loga=False, thres=2 )
    7099        sizey = detsize( yy, res = 10, loga=False, thres=2 )
     100        #sizex = detsize( xx, res = 10, loga=False, thres=3 )
     101        #sizey = detsize( yy, res = 10, loga=False, thres=3 )
    71102        ###
     103        #if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex  ### un peu limite dans certains cas
     104        if (len(sizex) > len(sizey))     : sizey = sizex        ### plus fidele mais petit souci lorsque PBC
     105        elif (len(sizex) == len(sizey))  :
     106          if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex
     107          else                                  : sizex = sizey
     108        else                            : sizex = sizey
    72109        allsizesx = np.append(allsizesx,sizex)
    73110        allsizesy = np.append(allsizesy,sizey)
     111        print i, ' on ', shape[0], ' caught ', len(sizex), ' vortices ', sizex
    74112        ########
    75    
     113 
    76114    allsizesx.sort()
    77115    allsizesy.sort()
    … …  
    79117    return allsizesx, allsizesy
    80118
    81 def writeascii ( tab, filename ):
    82     mydata = tab
    83     myfile = open(filename, 'w')
    84     for line in mydata:
    85         myfile.write(str(line) + '\n')
    86     myfile.close()
    87     return
     119#########################################################################
     120#########################################################################
    88121
    89122import matplotlib.pyplot as plt
    … …  
    91124import numpy as np
    92125import matplotlib.mlab as mlab
     126import mymath as mym
    93127
    94128save = True
    95 #save = False
     129save = False
     130
    96131if save:
    97132    allsizesx, allsizesy = getsize("psfc.nc")
    98     writeascii(allsizesx,'allsizex.txt')
    99     writeascii(allsizesy,'allsizey.txt')
    100    
     133    ### sauvegarde texte pour inspection
     134    mym.writeascii(allsizesx,'allsizex.txt')
     135    mym.writeascii(allsizesy,'allsizey.txt')
     136    ### sauvegarde binaire pour utilisation python
    101137    myfile = open('allsizex.bin', 'wb')
    102138    pickle.dump(allsizesx, myfile)
    103139    myfile.close()
    104    
    105140    myfile = open('allsizey.bin', 'wb')
    106141    pickle.dump(allsizesy, myfile)
    107142    myfile.close()
    108143
    109 
     144### load files
    110145myfile = open('allsizex.bin', 'r')
    111146allsizesx = pickle.load(myfile)
    … …  
    113148allsizesy = pickle.load(myfile)
    114149
     150### append sizes
    115151plothist = np.append(allsizesx,allsizesy)
     152plothist = allsizesx
     153#plothist = allsizesy
    116154plothist.sort()
    117155
    118 nbins=100
     156### make bins
     157nbins = 100
    119158zebins = [2.0]
     159#nbins = 8
     160#zebins = [19.0]
     161#nbins = 15
     162#zebins = [11.] ##12 non mais donne un peu la meme chose
     163#zebins = [20.] ##20 non car trop pres du premier
     164nbins = 100
     165zebins = [2./np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
     166nbins = 200
     167zebins = [0.2/np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
     168
    120169for i in range(0,nbins):  zebins.append(zebins[i]*np.sqrt(2))
    121170zebins = np.array(zebins)
    122 print zebins
    123 
     171### select reasonable bins for DD
    124172zebins = zebins [ zebins > 10. ]
    125173zebins = zebins [ zebins < 1000. ]
    126 print zebins
     174#print 'bins ',zebins
     175print 'corrected bins ',zebins
     176print 'max value ',mym.max(plothist),mym.max(allsizesx),mym.max(allsizesy)
     177
     178#plothist = [20.,30.,40.,50.,70.,100.,130.,190.,260.,370.,520.]  ## pour calibrer
    127179
    128180#### HISTOGRAM
    … …  
    130182             log=True,\
    131183             bins=zebins,\
    132              #cumulative=-1,\
     184#             cumulative=-1,\
    133185        )
    134186plt.xscale('log')
    135 
    136187plt.show()
    137188

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py

    r261 r310  
    1414        return u'\u00b0'
    1515
     16def writeascii ( tab, filename ):
     17    mydata = tab
     18    myfile = open(filename, 'w')
     19    for line in mydata:
     20        myfile.write(str(line) + '\n')
     21    myfile.close()
     22    return
     23

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py

    r296 r310  
    450450    from mymath import max,min,mean
    451451    ### might be convenient to add the missing value in arguments
    452     what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7)
     452    #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
     453    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7)
     454    else:          what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7)
    453455    print "new min ", min(what_I_plot)
    454456    what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
    455     what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
     457    what_I_plot[ what_I_plot > zevmax ] = zevmax
    456458    print "new max ", max(what_I_plot)
     459   
    457460    return what_I_plot
    458461
    … …  
    508511             "USTM":         "YlOrRd",\
    509512             "HFX":          "RdYlBu",\
    510              "ICETOT":       "YlGnBu",\
     513             "ICETOT":       "YlGnBu_r",\
    511514             "MTOT":         "PuBu",\
    512515             "TAU_ICE":      "Blues",\

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F

    r309 r310  
    497497      ztoploc = 10      ! target pseudo-altitude of local storm (km)
    498498      radloc = 0.5      ! radius of dust storm (degree)
    499       lonloc = 80       ! center longitude of storm (deg)
    500       latloc = -25      ! center latitude of storm (deg)
     499      lonloc = 25.      ! center longitude of storm (deg)
     500      latloc = -3.      ! center latitude of storm (deg)
    501501
    502502      DO ig=1,ngrid

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/initracer.F

    r308 r310  
    1 link ../../../mars_lmd_new/libf/phymars/initracer.F
     1link STORM_JULIEN_LAST/initracer.F
Note: See TracChangeset for help on using the changeset viewer.