Changeset 186


Ignore:
Timestamp:
Jul 3, 2011, 4:27:24 AM (14 years ago)
Author:
aslmd
Message:

MESOSCALE: python post-processing. wrapper with my Fortran interpolator (small modifications done to api.F90). added the corresponding options to winds.py which is now a pretty complete script

Location:
trunk/MESOSCALE
Files:
2 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/api.F90

    r112 r186  
    3535      INCLUDE 'netcdf.inc'
    3636
     37      !
     38      ! VARIABLES
     39      !
     40      CHARACTER (LEN=500)                                :: path_to_input
     41      CHARACTER (LEN=500)                                :: path_to_output
     42      CHARACTER (LEN=500)                                :: input_name
     43      CHARACTER (LEN=500)                                :: output_name
     44      CHARACTER (LEN=20)                                 :: process
     45      CHARACTER (LEN=2000)                               :: fields
     46      REAL, DIMENSION(299)                               :: interp_levels
     47      INTEGER                                            :: interp_method=1
     48      INTEGER                                            :: extrapolate=0
     49      LOGICAL                                            :: debug=.FALSE.
     50      LOGICAL                                            :: unstagger_grid=.FALSE.
     51      LOGICAL                                            :: bit64=.FALSE.
     52      LOGICAL                                            :: oldvar=.TRUE.             
     53
     54      INTEGER                                            :: funit,ios
     55      LOGICAL                                            :: is_used
     56
     57      !
     58      ! NAMELISTS
     59      !
     60      NAMELIST /io/ path_to_input, input_name, path_to_output, output_name, &
     61                    process, fields, debug, bit64, oldvar
     62      NAMELIST /interp_in/ interp_levels, interp_method, extrapolate, unstagger_grid
     63 
     64      !
     65      ! DEFAULT VALUES for VARIABLES
     66      !
     67      path_to_input   = './'
     68      path_to_output  = './'
     69      output_name     = ' '
     70      interp_levels   = -99999.
     71      process         = 'all'
     72     
     73
     74      !
     75      ! READ NAMELIST
     76      !
     77        DO funit=10,100
     78           INQUIRE(unit=funit, opened=is_used)
     79           IF (.not. is_used) EXIT
     80        END DO
     81        OPEN(funit,file='namelist.api',status='old',form='formatted',iostat=ios)
     82        IF ( ios /= 0 ) STOP "ERROR opening namelist.api"
     83        READ(funit,io)
     84        READ(funit,interp_in)
     85        CLOSE(funit)
     86
     87      !!! MAIN CALL
     88      CALL api_main ( path_to_input, input_name, path_to_output, output_name, &
     89                                 process, fields, debug, bit64, oldvar, &
     90                                 interp_levels, interp_method, extrapolate, unstagger_grid, -99999. )
     91
     92 END PROGRAM api
     93
     94 SUBROUTINE api_main ( path_to_input, input_name, path_to_output, output_name, &
     95                       process, fields, debug, bit64, oldvar, &
     96                       interp_levels, interp_method, extrapolate, unstagger_grid, onelevel )
     97
     98      IMPLICIT NONE
     99      INCLUDE 'netcdf.inc'
     100
    37101      !!
    38102      !! EARTH CONSTANTS
    … …  
    62126      CHARACTER,         ALLOCATABLE, DIMENSION(:,:,:,:) :: text
    63127      CHARACTER (LEN=31),ALLOCATABLE, DIMENSION(:)       :: dnamei, dnamej
    64       CHARACTER(LEN=250),ALLOCATABLE, DIMENSION(:)       :: input_file_names
    65       CHARACTER(LEN=250),ALLOCATABLE, DIMENSION(:)       :: output_file_names
     128      CHARACTER(LEN=500),ALLOCATABLE, DIMENSION(:)       :: input_file_names
     129      CHARACTER(LEN=500),ALLOCATABLE, DIMENSION(:)       :: output_file_names
    66130      DOUBLE PRECISION,  ALLOCATABLE, DIMENSION(:,:,:,:) :: ddata1, ddata2
    67131      REAL,              ALLOCATABLE, DIMENSION(:,:,:,:) :: data1, data2, data3
    … …  
    79143      INTEGER,                        DIMENSION(4)       :: dims_in, dims_out
    80144      INTEGER,                        DIMENSION(6)       :: ishape, jshape
    81       CHARACTER (LEN=80)                                 :: cval
     145      CHARACTER (LEN=500)                                :: cval !80
    82146      CHARACTER (LEN=31)                                 :: cname, test_dim_name
    83       CHARACTER (LEN=80)                                 :: input_file, output_file, att_text
    84       CHARACTER (LEN=250)                                :: path_to_input
    85       CHARACTER (LEN=250)                                :: path_to_output
    86       CHARACTER (LEN=250)                                :: input_name
    87       CHARACTER (LEN=250)                                :: output_name, tmp_name
     147      CHARACTER (LEN=500)                                :: input_file, output_file, att_text !80
     148      CHARACTER (LEN=500)                                :: path_to_input
     149      CHARACTER (LEN=500)                                :: path_to_output
     150      CHARACTER (LEN=500)                                :: input_name
     151      CHARACTER (LEN=500)                                :: output_name, tmp_name
    88152      CHARACTER (LEN=10)                                 :: option
    89153      CHARACTER (LEN=132)                                :: command
    90154      CHARACTER (LEN=20)                                 :: process, dummy
    91155      CHARACTER (LEN=2000)                               :: fields, process_these_fields
    92       REAL, DIMENSION(299)                               :: interp_levels
     156      REAL, DIMENSION(299)                               :: interp_levels 
    93157      REAL                                               :: rval
    94158      REAL                                               :: MISSING=1.e36
    … …  
    96160      INTEGER                                            :: map_proj
    97161      INTEGER                                            :: LINLOG = 1
    98       INTEGER                                            :: interp_method=1
    99       INTEGER                                            :: extrapolate=0
     162      INTEGER                                            :: interp_method!=1
     163      INTEGER                                            :: extrapolate!=0
    100164      INTEGER                                            :: ncid, mcid, rcode
    101165      INTEGER                                            :: idm, ndims, nvars, natt, ngatts
    … …  
    111175      INTEGER                                            :: kk
    112176      LOGICAL                                            :: is_used
    113       LOGICAL                                            :: debug=.FALSE.
     177      LOGICAL                                            :: debug!=.FALSE.
    114178      LOGICAL                                            :: interpolate=.FALSE.
    115       LOGICAL                                            :: unstagger_grid=.FALSE.
     179      LOGICAL                                            :: unstagger_grid!=.FALSE.
    116180      LOGICAL                                            :: fix_meta_stag=.FALSE.
    117       LOGICAL                                            :: bit64=.FALSE.
     181      LOGICAL                                            :: bit64!=.FALSE.
    118182      LOGICAL                                            :: first=.TRUE.
    119       LOGICAL                                            :: oldvar=.FALSE.             
    120 
    121       !
    122       ! NAMELISTS
    123       !
    124       NAMELIST /io/ path_to_input, input_name, path_to_output, output_name, &
    125                     process, fields, debug, bit64, oldvar
    126       NAMELIST /interp_in/ interp_levels, interp_method, extrapolate, unstagger_grid
    127 
    128       !
    129       ! DEFAULT VALUES for VARIABLES
    130       !
    131       path_to_input   = './'
    132       path_to_output  = './'
    133       output_name     = ' '
    134       interp_levels   = -99999.
    135       process         = 'all'
    136 
    137       !
    138       ! READ NAMELIST
    139       !
    140         DO funit=10,100
    141            INQUIRE(unit=funit, opened=is_used)
    142            IF (.not. is_used) EXIT
    143         END DO
    144         OPEN(funit,file='namelist.api',status='old',form='formatted',iostat=ios)
    145         IF ( ios /= 0 ) STOP "ERROR opening namelist.api"
    146         READ(funit,io)
    147         READ(funit,interp_in)
    148         CLOSE(funit)
     183      LOGICAL                                            :: oldvar!=.FALSE.                   
     184
     185      REAL :: onelevel
     186      if ( onelevel .ne. -99999. ) then
     187        interp_levels(1) = onelevel
     188        interp_levels(2:) = -99999.
     189      endif
    149190
    150191      !
    … …  
    152193      !
    153194        lent = len_trim(path_to_input)
    154         IF ( path_to_input(lent:lent) /= "/" ) THEN
    155            path_to_input = TRIM(path_to_input)//"/"
    156         ENDIF
    157         lent = len_trim(path_to_output)
    158         IF ( path_to_output(lent:lent) /= "/" ) THEN
    159            path_to_output = TRIM(path_to_output)//"/"
    160         ENDIF
     195         !IF ( path_to_input(lent:lent) /= "/" ) THEN
     196         !   path_to_input = TRIM(path_to_input)//"/"
     197         !ENDIF
     198         !lent = len_trim(path_to_output)
     199         !IF ( path_to_output(lent:lent) /= "/" ) THEN
     200         !   path_to_output = TRIM(path_to_output)//"/"
     201         !ENDIF
    161202        input_name = TRIM(path_to_input)//TRIM(input_name)
    162203        !
    … …  
    526567        write(6,*) "  Data will be output on unstaggered grid "
    527568        do kk = 1, times_in_file
    528          IF ( DEBUG ) print *, kk
     569         !IF ( DEBUG ) print *, kk
    529570         IF (oldvar) THEN
    530571           interm1(1:iweg-1,:,:) = ( u(1:iweg-1,:,:,kk) + u(2:iweg,:,:,kk) ) * .5
    … …  
    673714                ENDDO
    674715                   deallocate (data1)
    675                    PRINT *, pres_field(10,10,:,1)
     716                   !PRINT *, pres_field(10,10,:,1)
    676717           ENDIF
    677718 
    … …  
    949990                                      DO ii = 1,4
    950991                                        dims_out(ii) = dvalj(jshape(ii))
    951                                         print *, dims_out(ii)
     992                                        !print *, dims_out(ii)
    952993                                      ENDDO
    953994                                      !!! NB: what follows is useful because we'd like diagnostics for each history timestep
    … …  
    11341175      write(6,*) "##########################################"
    11351176
    1136  END PROGRAM api
     1177END SUBROUTINE
     1178! END PROGRAM api
    11371179!---------------------------------------------------------------------
    11381180!---------------------------------------------------------------------
    … …  
    15731615   real, parameter :: RAD_PER_DEG = PI/180.
    15741616
    1575 !  print *, 'map ', map_proj
    15761617  IF ( map_proj .ge. 3 ) THEN                         ! No need to rotate
    15771618    !PRINT *, 'NO NEED TO ROTATE !!!! equivalent to output U,V with unstagger_grid'
    … …  
    15791620    VVVmet(:,:,:) = VVV
    15801621  ELSE
    1581   !END IF
    15821622
    15831623  cone = 1.                                          !  PS
    … …  
    16061646  END DO
    16071647
    1608 !print *, longi(10,10)
    1609 !print *, lati(10,10)
    1610 
    16111648
    16121649  DO i = 1, west_east_dim
    … …  
    16241661    UUUmet(:,:,k) = VVV(:,:,k)*sin(alpha) + UUU(:,:,k)*cos(alpha)
    16251662    VVVmet(:,:,k) = VVV(:,:,k)*cos(alpha) - UUU(:,:,k)*sin(alpha)
    1626 !print *,UUU(10,10,k), UUUmet(10,10,k)
    1627 !print *,VVV(10,10,k), VVVmet(10,10,k)
    16281663  END DO
    16291664  END IF

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/comp_api

    r114 r186  
    55#
    66
    7 pgf90 api.F90 -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include -Mfree -o api
     7#pgf90 api.F90 -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include -Mfree -o api
    88
    99
    10 #g95 api.F90 -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include -Wall -Wno=112,141,137,155 -fno-second-underscore -ffree-form -o api
     10g95 api.F90 -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include -Wall -Wno=112,141,137,155 -fno-second-underscore -ffree-form -o api
    1111#g95 api.F90 -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include -Wall -Wno=112,141,137,155 -ffree-form -o api
    1212#pgf90 -mcmodel=medium -Mlarge_arrays api.F90 -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include -Mfree -o p_interp

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/namelist.api

    r19 r186  
    11&io
    2 ! path_to_input  = './dossier/ALT_POLAR_top0.1_sponge_smooth_topo/'  !! ne doit pas etre trop long
    3 ! path_to_output = './smooth'
    4  path_to_input  = './dossier/ALT_POLAR_top0.1_sponge/'
     2 path_to_input  = './'
    53 path_to_output = './'
    6 ! input_name     = 'wrfout_d01_2024*'
    7  input_name     = 'wrfout_d01_2024-07-03_06:00:00'
    8 ! path_to_input = './dossier_test/'
    9 ! input_name     = 'wrfout_d01_2024-01-05_02:00:00'
    10 ! input_name     = 'nca_wrfout_d01_2024-01-05_02:00:00'
     4 input_name     = 'wrfout_d01_9999-09-09_09:00:00'
    115 process        = 'list'    !! list fields required in "fields" (available tk, tpot, GHT)     
    12 ! fields         = 'W,SWDOWNZ,TAU_DUST,TSURF,XLONG,XLAT,HGT'
    13 ! fields         = 'W'
    14 ! fields         = 'U,V'  !! beware, not only tk ! only one 3D field at least
    156 fields         = 'tk,W,uvmet' 
    16  debug          = .TRUE.
    177/
    18  process = 'all'  !! process all fields in file
    19  debug = .TRUE.
    20  bit64 = .TRUE. ! bit64 = .TRUE. !! ne pas utiliser ?
    218
    229&interp_in
    23 ! interp_levels = 1000.,987.5,975.,962.5,950.,937.5,925.,
    24 !                 900.,875.,850.,825.,800.,750.,700.,650., 
    25 !                 600.,550.,500.,450.,400.,350.,300.,250.,
    26 !                 225.,200.,175.,150.,137.5,125.,112.5,100.,
    27 !                 87.5,75.,62.5,50.,37.5,25.,12.5,
    28 !interp_levels = 5., 4.9, 4.8, 4.7, 4.6, 4.5, 4.4, 4.3, 4.2, 4.1, 4.0,
    29 !                    3.9, 3.8, 3.7, 3.6, 3.5, 3.4, 3.3, 3.2, 3.1, 3.0,
    30 !                    2.9, 2.8, 2.7, 2.6, 2.5, 2.4, 2.3, 2.2, 2.1, 2.0,
    31 !                    1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1, 1.0,
    32 !                    0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1,
    33 ! interp_levels = 5., 3., 1., 0.1
    34 ! interp_method = 2,
    35  interp_method = 3                 !! INTERPOLATION: ALTITUDE ABOVE MOLA AREOID (km)
    36 ! interp_method = 4
    37 ! interp_levels = 20.
    38 ! interp_levels = -09., -08., -07., -06., -05., -04., -03., -02., -01.,  00.,
    39 !                  01.,  02.,  03.,  04.,  05.,  06.,  07.,  08.,  09.,
    40 ! interp_levels =  01.,  02.,  03.,  04.,  05.,  06.,  07.,  08.,  09.,
    41 !                  10.,  11.,  12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,
    42 !                  20.,  21.,  22.,  23.,  24.,  25.,  26.,  27.,  28.,  29.,
    43 !!                  30.,  31.,  32.,  33.,  34.,  35.,  36.,  37.,  38.,  39.,
    44 !!                  40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,  48.,  49.
    45 ! interp_levels =               01.,  01.5,  02.,  02.5,  03.,  03.5,  04.,  04.5,  05.,  05.5,  06.,  06.5,  07.,  07.5,  08.,  08.5,  09.,  09.5
    46 !                  10.,  10.5,  11.,  11.5,  12.,  12.5,  13.,  13.5,  14.,  14.5,  15.,  15.5,  16.,  16.5,  17.,  17.5,  18.,  18.5,  19.,  19.5
    47 !                  20.,  20.5,  21.,  21.5,  22.,  22.5,  23.,  23.5,  24.,  24.5,  25.,  25.5,  26.,  26.5,  27.,  27.5,  28.,  28.5,  29.,  29.5
    48  interp_levels = 20.0
    49 ! interp_levels = 01.00,01.25,01.50,01.75,02.00,02.25,02.50,02.75,03.00,03.25,03.50,03.75,04.00,04.25,04.50,04.75,05.00,
    50 !                       05.25,05.50,05.75,06.00,06.25,06.50,06.75,07.00,07.25,07.50,07.75,08.00,08.25,08.50,08.75,09.00,
    51 !                       09.25,09.50,09.75,10.00,10.25,10.50,10.75,11.00,11.25,11.50,11.75,12.00,12.25,12.50,12.75,13.00,
    52 !                       13.25,13.50,13.75,14.00,14.25,14.50,14.75,15.00,15.25,15.50,15.75,16.00,16.25,16.50,16.75,17.00
    53 ! interp_method = 2,
    54 ! interp_levels = 1.,
     10 interp_method = 4
     11 interp_levels = 0.050
    5512/
    56  extrapolate = 1,
    57  interp_method = 2,
    58  unstagger_grid = .TRUE.  !! pb !! non c OK
    59 
    60  extrapolate = 0    ;; set values below ground and above model top to missing (default)
    61  extrapolate = 1    ;; extrapolate below ground, and set above model top to model top values
    62  interp_method = 1  ;; linear in p interpolation (default)
    63  interp_method = 2  ;; linear in log p interpolation

  • trunk/MESOSCALE/PLOT/PYTHON/mylib/myplot.py

    r184 r186  
    3535        return wlon,wlat
    3636
     37def api_onelevel (  path_to_input   = None, \
     38                    input_name      = 'wrfout_d0?_????-??-??_??:00:00', \
     39                    path_to_output  = None, \
     40                    output_name     = 'output.nc', \
     41                    process         = 'list', \
     42                    fields          = 'tk,W,uvmet,HGT', \
     43                    debug           = False, \
     44                    bit64           = False, \
     45                    oldvar          = True, \
     46                    interp_method   = 4, \
     47                    extrapolate     = 0, \
     48                    unstagger_grid  = False, \
     49                    onelevel        = 0.020  ):
     50    import api
     51    import numpy as np
     52    if not path_to_input:   path_to_input  = './'
     53    if not path_to_output:  path_to_output = path_to_input
     54    api.api_main ( path_to_input, input_name, path_to_output, output_name, \
     55                   process, fields, debug, bit64, oldvar, np.arange (299), \
     56                   interp_method, extrapolate, unstagger_grid, onelevel )
     57    return
     58
    3759def getproj (nc):
    3860    map_proj = getattr(nc, 'MAP_PROJ')
    … …  
    5173        print "mercator projection"
    5274        proj="merc"
     75    else:
     76        proj="merc"
    5377    return proj   
    5478
    … …  
    7094    import  matplotlib.pyplot as plt
    7195    res = int(res)
    72     if folder != '':      name = folder+'/'+filename+str(res)+".png"
    73     else:                 name = filename+str(res)+".png"
     96    name = filename+"_"+str(res)+".png"
     97    if folder != '':      name = folder+'/'+name
    7498    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
    7599    if disp:              display(name)         

  • trunk/MESOSCALE/PLOT/PYTHON/scripts/winds.py

    r184 r186  
    2424    import numpy as np
    2525
    26     #############################
    27     ### Lower a bit the font size
    28     rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
    29 
    3026    ######################
    3127    ### Load NETCDF object
    … …  
    5046                                    ## pb avec les autres (de trace derriere la sphere ?)
    5147
    52     ###################################################################
    53     ### For mesoscale results plot the underlying topography by default
     48    ####################################################################
     49    #### For mesoscale results plot the underlying topography by default
    5450    if typefile in ['mesoapi','meso']:
    55         if var == None: var = 'HGT'
     51        if var == None: back="mola" #var = 'HGT'
    5652
    5753    ####################################################
    … …  
    8379        [u,v] = getwinds(nc,charu='U',charv='V')
    8480        metwind = False ## geometrical (wrt grid)
     81        print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
    8582
    8683    ####################################################
    … …  
    9592            elif dimension == 3:   field = nc.variables[var][:,:,:]
    9693            elif dimension == 4:   field = nc.variables[var][:,nvert,:,:] 
     94    nt = len(u[:,0,0,0])
     95
     96    #########################################
     97    ### Name for title and graphics save file
     98    basename = "WINDS"
     99    if var:
     100        basename = basename + "_" + var
     101    basename = basename + "_z" + str(nvert)
    97102
    98103    ##################################
    99104    ### Open a figure and set subplots
    100105    fig = figure()
    101     if   numplot == 4:
    102         sub = 221
    103         fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
    104     elif numplot == 2:
    105         sub = 121
    106         fig.subplots_adjust(wspace = 0.3)
    107     elif numplot == 3:
    108         sub = 131
    109         fig.subplots_adjust(wspace = 0.3)
    110     elif numplot == 6:
    111         sub = 231
    112         fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
    113     elif numplot == 8:
    114         sub = 331 #241
    115         fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
    116     else:
    117         print "supported: 1,2,3,4,6,8"
    118         exit()
    119 
    120     #####################
    121     ### Prepare time loop
    122     nt = len(u[:,0,0,0])
    123     if nt <= numplot or numplot == 1: 
    124         print "I am plotting only one map ",nt,numplot
    125         tabrange = [0]
    126     else:                         
    127         tabrange = range(0,nt-1,int(nt/numplot))
     106    if   numplot > 0:   
     107        if   numplot == 4:
     108            sub = 221
     109            fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
     110            rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
     111        elif numplot == 2:
     112            sub = 121
     113            fig.subplots_adjust(wspace = 0.3)
     114            rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
     115        elif numplot == 3:
     116            sub = 131
     117            fig.subplots_adjust(wspace = 0.3)
     118            rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
     119        elif numplot == 6:
     120            sub = 231
     121            fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
     122            rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
     123        elif numplot == 8:
     124            sub = 331 #241
     125            fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
     126            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
     127        elif numplot == 9:
     128            sub = 331
     129            fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
     130            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
     131        elif numplot == 1:
     132            pass
     133        else:
     134            print "supported: 1,2,3,4,6,8"
     135            exit()
     136        ### Prepare time loop
     137        if nt <= numplot or numplot == 1: 
     138            tabrange = [0]
     139            numplot = 1
     140        else:                         
     141            tabrange = range(0,nt,int(nt/numplot))  #nt-1
     142            tabrange = tabrange[0:numplot]
     143    else:
     144        tabrange = range(0,nt,1)
     145        sub = 99999
     146    print tabrange
    128147
    129148    #################################
    130149    ### Time loop for plotting device
     150    found_lct = False
    131151    for i in tabrange:
    132152
    133        print i
    134 
    135153       ### General plot settings
    136        if tabrange != [0]: subplot(sub)
    137        zetitle = "WINDS" + "_"
     154       if numplot > 1:
     155           subplot(sub)
     156           found_lct = True
     157       elif numplot == 1:
     158           found_lct = True
     159        ### If only one local time is requested (numplot < 0)
     160       elif numplot <= 0:
     161           #print (ltst+i)%24, numplot, (ltst+i)%24+numplot
     162           if (ltst+i)%24 + numplot != 0:   continue
     163           else:                            found_lct = True
    138164
    139165       ### Map projection
    … …  
    143169       #### Contour plot
    144170       if var:
    145            zetitle = zetitle + var + "_"
    146171           if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(field[i,:,:])
    147172           elif typefile in ['gcm']:             
    … …  
    158183           key = False
    159184       if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
    160        else:        print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
    161185       vectorfield(vecx, vecy,\
    162                       x, y, stride=stride, csmooth=2,\
     186                      x, y, stride=stride, csmooth=stride,\
    163187                      scale=15., factor=200., color='k', key=key)
     188                                                   ## or csmooth=2
    164189       
    165190       ### Next subplot
    166        zetitle = zetitle + "LT"+str((ltst+i)%24)
    167        ptitle(zetitle)
     191       ptitle( basename + "_LT"+str((ltst+i)%24) )
    168192       sub += 1
    169193
    170194    ##########################################################################
    171195    ### Save the figure in a file in the data folder or an user-defined folder
    172     if not target:   zeplot = namefile+zetitle
    173     else:            zeplot = target+"/"+zetitle
    174     makeplotpng(zeplot,pad_inches_value=0.35)   
     196    if not target:    zeplot = namefile+"_"+basename
     197    else:             zeplot = target+"/"+basename
     198    if numplot <= 0:  zeplot = zeplot + "_LT"+str(abs(numplot))
     199    if found_lct:     makeplotpng(zeplot,pad_inches_value=0.35)   
     200    else:             print "Local time not found"
    175201
    176202
    … …  
    196222    import sys
    197223    from optparse import OptionParser    ### to be replaced by argparse
     224    from api_wrapper import api_onelevel
    198225    parser = OptionParser()
    199226    parser.add_option('-f', action='store', dest='namefile',    type="string",  default=None,  help='name of WRF file [NEEDED]')
    200     parser.add_option('-l', action='store', dest='nvert',       type="int",     default=0,     help='subscript for vertical level')
     227    parser.add_option('-l', action='store', dest='nvert',       type="float",   default=0,     help='vertical level (def=0)')
    201228    parser.add_option('-p', action='store', dest='proj',        type="string",  default=None,  help='projection')
    202229    parser.add_option('-b', action='store', dest='back',        type="string",  default=None,  help='background')
    203230    parser.add_option('-t', action='store', dest='target',      type="string",  default=None,  help='destination folder')
    204     parser.add_option('-s', action='store', dest='stride',      type="int",     default=3,     help='stride vectors')
     231    parser.add_option('-s', action='store', dest='stride',      type="int",     default=3,     help='stride vectors (def=3)')
    205232    parser.add_option('-v', action='store', dest='var',         type="string",  default=None,  help='variable contoured')
    206     parser.add_option('-n', action='store', dest='numplot',     type="int",     default=4,     help='number of plots')
     233    parser.add_option('-n', action='store', dest='numplot',     type="int",     default=4,     help='number of plots (def=1)(if <0: 1 plot of LT -*numplot*)')
     234    parser.add_option('-i', action='store', dest='interp',      type="int",     default=None,  help='interpolation method (done at level *nvert* km)')
    207235    (opt,args) = parser.parse_args()
    208236    if opt.namefile is None:
    … …  
    210238        exit()
    211239    print "Options:", opt
    212     winds (opt.namefile,opt.nvert,proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot)
    213 
    214 #    if typefile in ['gcm']:
    215 #        if var == 'HGT':    var = 'phisinit'  ## default choice for GCM
    216 
    217 
     240
     241    zefile = opt.namefile   
     242    zelevel = opt.nvert   
     243    if opt.nvert is 0 and opt.interp:   zelevel = 0.020
     244    if opt.interp is not None:
     245        if   opt.var is None    :  zefields = 'uvmet'
     246        else                    :  zefields = 'uvmet,'+opt.var
     247        zefile = api_onelevel (  path_to_input   = '', \
     248                                 input_name      = zefile, \
     249                                 path_to_output  = opt.target, \
     250                                 fields          = zefields, \
     251                                 interp_method   = opt.interp, \
     252                                 onelevel        = zelevel )
     253        zelevel = 0
     254
     255    winds (zefile,int(zelevel),proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot)
Note: See TracChangeset for help on using the changeset viewer.