source: trunk/MESOSCALE_DEV/PLOT/RESERVE/out_geo.pro @ 603

pro out_geo, $
        field1=charvar, $
        field2=charvar2, $
        domain=domain, $
        path=path

;---------------------------------------------------------------;
;                                                               ;
; Use:                                                          ;
;       out_geo, 'HGT_M'        (or) out_geo, 'TOPO'            ;
;       out_geo, 'THERMAL_INERTIA' (or) out_geo, 'TI'           ;
;       out_geo, 'ALBEDO_GCM' (or) out_geo, 'ALBEDO'            ;
;       out_geo, 'GHT'                                          ;
;                                                               ;
;       -- nested domain                                        ;
;       out_geo, 'HGT_M', domain='d02'                          ;
;                                                               ;
; Options:                                                      ;
;               - domain: default is 'd01',                     ;
;                       but can be 'd02', 'd03', etc ...        ;
;               - path: where the geo_em files are              ; 
;                       '/my_root/my_folder/'                   ;
;                                                               ;
; A. Spiga, April 2007 - July 2007 - September 2007             ;
;                                                               ;
;---------------------------------------------------------------;




;----------------------
; set parameters
;----------------------

colors=32
format='(F6.2)'
minfield_init=0.
maxfield_init=0.

default_path='./'
if (n_elements(path) eq 0) then path=default_path
if (n_elements(domain) eq 0) then domain='d01'
file=path+'/geo_em.'+domain+'.nc'



;----------------------
; equivalent inputs
;----------------------
case charvar of
        'HGT_M':
        'TOPO': charvar='HGT_M'
        'THERMAL_INERTIA':
        'TI': charvar='THERMAL_INERTIA'
        'ALBEDO_GCM':
        'ALBEDO': charvar='ALBEDO_GCM'
        'GHT':
        else:
endcase 
if (n_elements(charvar2) ne 0) then begin
case charvar2 of
        'HGT_M':
        'TOPO': charvar2='HGT_M'
        'THERMAL_INERTIA':
        'TI': charvar2='THERMAL_INERTIA'
        'ALBEDO_GCM':
        'ALBEDO': charvar2='ALBEDO_GCM'
        'GHT':
        else:
endcase
endif                                                   


;----------------------
; set graphics
;----------------------

set_plot, 'ps'
device, file='geo_em.'+domain+'_'+charvar+'.ps', /color
;set_plot, 'x'

case charvar of
        'HGT_M': pal=33 ;16  ;4
        'THERMAL_INERTIA': pal=3
        'ALBEDO_GCM': pal=0
        'GHT': pal=33
        else:
endcase
                                


!p.charthick = 2.0
!p.thick = 3.0
!x.thick = 2.0
!y.thick = 2.0

case charvar of
        'HGT_M': title='Altimetry (km)'
        'THERMAL_INERTIA': title='Thermal inertia (!NJ!N.m!U-2!N.s!U-1/2!N.K!U-1!N)'
        'ALBEDO_GCM': title='Albedo LW (%)'
        'GHT': title='Geopotential height (m)'
        else: title=''
endcase


;-------------------------------
; open file and read variables
;-------------------------------
cdfid = ncdf_open(file)

; field
print, charvar
varid=ncdf_varid(cdfid,charvar)
        ncdf_varget, cdfid, varid, var
; lon
varid=ncdf_varid(cdfid,'XLONG_M')
        ncdf_varget, cdfid, varid, lon
        lon=reform(lon[*,0])
; lat
varid=ncdf_varid(cdfid,'XLAT_M')
        ncdf_varget, cdfid, varid, lat
        lat=reform(lat[0,*])

if (n_elements(charvar2) ne 0) then begin
        print, charvar2
        varid=ncdf_varid(cdfid,charvar2)
        ncdf_varget, cdfid, varid, overcontour
endif else begin
        overcontour=0.
endelse 
        

;-------------------------------
; calculation on variables
;-------------------------------
case charvar of
        'HGT_M': var=var/1000
        'THERMAL_INERTIA': 
        'ALBEDO_GCM': var=var*100
        'GHT': 
        else:
endcase
                                
        
        
;---------
; plot !
;---------

minfield=min(var)
maxfield=max(var)
print, 'mM', minfield, maxfield

what_I_plot=var


map_set, -90, 0, /azimuthal, limit=[-90,-180,-60,180]
map_latlon, $
        what_I_plot, $                          ; 2D field
        lon, $                                  ; 1D latitude
        lat                                  ; 1D longitude
 

map_latlon, $
        what_I_plot, $                          ; 2D field
        lon, $                                  ; 1D latitude
        lat, $                                  ; 1D longitude
        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
        overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
;        overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
;        overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
        ct=pal, $                               ; color table (33-rainbow is default)
        colors=colors, $                        ; number of colors/levels (32 is default)
        title=title, $                          ; title of the plot ('' is default)
        format=format                           ; format of colorbar annotations ('(F6.2)' is default)
                                                                                                

device, /close

end