source: trunk/mesoscale/PLOT/RESERVE/obsolete/out_wrf.pro @ 106

pro out_wrf, $
                plot=plot_user, $               ; 'map', 'profile', 'hovmuller', 'meridional', 'zonal' (default is 'map')
        ;----------------
        ; input fields
        ;----------------
                field1=field_user, $            ; field to contour: 'U', 'V', 'tk', ... etc ...
                field2=fieldc_user, $           ; field to overplot (contour lines, optional)
        ;------------------------------------------------------
        ; 'map' options (no effect when plot is 'profile')
        ;------------------------------------------------------
                topo=topo, $                    ; overplot topography contours (!! must load HGT !! override field2)
                winds=winds, $                  ; option to overplot winds (e.g. =['U','V'])
                level=level, $                  ; vertical level subscript (starts at 1 ...)(if -1: all levels)
        ;------------------------------------------------------
        ; 'profile' options (no effect when plot is 'map')
        ;------------------------------------------------------
                inprofile=inprofile, $          ; a profile to compare with ... (override field2)
                incolumn=incolumn, $            ; altitudes ... in case /= column
                ;; NB: inprofile and incolumn may be used to output profiles
                ;; ... reset with inprofile=0, incolumn=0
        ;----------------
        ; choose time
        ;----------------
                when=when, $                    ; time subscript (starts at 1 ...)(if -1: first time step)
        ;----------------
        ; choose latitude/longitude
        ;----------------
                nlat=nlat, $                    ; latitude subscript (starts at 1 ...)(if -1: all latitudes)(default is middle)
                nlon=nlon, $                    ; longitude subscript (starts at 1 ...)(if -1: all longitudes)(default is middle)
        ;----------------
        ; plot options
        ;----------------
;;NB: celles ci sont dorenavant indiquees dans le _user si besoin 
                range=range_user, $             ; range=[min,max]
                colors=colors, $                ; 16, 32, 64, ... etc
save_ps=file_ps, $
        ;-------------------------
        ; no useless data loading
        ;-------------------------
                save_data=save_data             ; to avoid loading data if already done once



;---------------------------------------------------------------------------
;
; A general routine to plot ARWpost outputs
;
; USE:
; - out_wrf     (will plot everything !)
; - out_wrf, field1='tk', range=[160,300]
;
; - out_wrf, field1='pressure', field2='tk'
; - out_wrf, field1=TSK, field2=HGT, winds=['U','V'] 
;
; - out_wrf, save_data=field
;   ... then next calls including 'save_data=field' won't include the data loading 
;
; - out_wrf, field1='tk', winds=['umet','vmet'], topo=1, save_data=yeah, level=1, when=2
;
; - out_wrf, field1='tk', plot='profile'
;
; - out_wrf, plot='profile', field1='tk', when=1, inprofile=inprofile, incolumn=incolumn
;   out_diagfi, plot='profile', field1='temp', when=112, inprofile=inprofile, incolumn=incolumn         
;
; - out_wrf, plot='meridional',field1='U'
;
; - out_wrf, plot='hovmuller',field1='TSK'
;
; - out_wrf, plot='zonal',field1='V', nlat=-1
; 
;---------------------------------------------------------------------------
; A. Spiga, August 2007, adding plot keyword September 2007
;---------------------------------------------------------------------------


    


;-------------------------------------
; default settings and user settings
;-------------------------------------
diff_scale=0  ;;1
charend=''
format=''
if (n_elements(plot_user) eq 0) then plot_user='map'
if (n_elements(field_user) eq 0) then field_user='all'
if (n_elements(level) eq 0) then level=-1
if (n_elements(when) eq 0) then when=-1
;;if (n_elements(nlon) eq 0) then nlon=-1
;;if (n_elements(nlat) eq 0) then nlat=-1
if (n_elements(colors) eq 0) then colors=32
;if (n_elements(minfield_init)*n_elements(maxfield_init) eq 0) then begin
if (n_elements(range_user) eq 0) then begin
        minfield_init = 0.
        maxfield_init = 0.
diff_scale=1
endif

if (n_elements(range_user) eq 2) then begin
        minfield_init=float(range_user(0))
        maxfield_init=float(range_user(1))
        print, 'setting user limits !', minfield_init, maxfield_init 
        diff_scale=1
endif

if (n_elements(topo) eq 0) then topo=0
if (topo eq 1) then fieldc_user='HGT'   ;case with topography

if (n_elements(fieldc_user) eq 0) then fieldc_user=''

;goto, squeeze
;------------------
; read data
;------------------
;
; read the file generated by ARWpost
;
;*********************************************
if (n_elements(save_data) eq 0) then begin
        read_dat, field, x, y, z, t, desc_field
        ;save_data=field
endif else begin
        read_dat, field, x, y, z, t, desc_field, nodata=1  ; read only the dimensions
        field=save_data
        save_data=0.    ;;free memory
endelse
;;*********************************************
;; opt: saved IDL files
;squeeze:
;restore, filename='OM_nonhydro.idl'
;field=output
;;nh=output
;;restore, filename='OM_hydro.idl'
;;h=output
;;field=abs(nh-h)
;;*********************************************

;
; dimensions
;
lon=x
west_east=n_elements(x)
if (n_elements(nlon) eq 0) then begin
        nlon=round(west_east/2)
endif 
if (nlon eq 0) then nlon=round(west_east/2)

lat=y
south_north=n_elements(y)
if (n_elements(nlat) eq 0) then begin
        nlat=round(south_north/2)
endif
if (nlat eq 0) then nlat=round(south_north/2)
        
bottom_top=n_elements(z)

time=n_elements(t)
if ((time eq 1) and (plot_user eq 'hovmuller')) then stop

nfields=n_elements(desc_field)
if (z(0) lt z(1)) then vert_coord='height'
if (z(0) gt z(1)) then vert_coord='pressure'
if ((z(0) eq 1.) and (z(bottom_top-1) eq (z(bottom_top-2)+1.))) then vert_coord='model_level'
;print, 'vert_coord is ... ', vert_coord

;
; find the field desired by the user
;
if (field_user ne 'all') then begin
        match=-1
        indice=-1
        repeat begin
                indice=indice+1
                if (indice eq nfields) $
                        then match=9999 $
                        else match=STRPOS(desc_field(indice),field_user)  
        endrep until match ne -1
        if (match eq 9999) then begin
                print, field_user, ' not found !'
                stop
        endif
;print, 'selected for plot:  ',desc_field(indice)
endif

;
; prepare loop for graphics
;
if (n_elements(indice) ne 0) then begin
        field_loop_start=indice
        field_loop_end=indice
endif else begin
        field_loop_start=0
        field_loop_end=nfields-1        
endelse

;
; overplot contour of another field ?
;
if (fieldc_user ne '') then begin
        match=-1
        indice=-1
        repeat begin
                indice=indice+1
                if (indice eq nfields) $
                        then match=9999 $
                        else match=STRPOS(desc_field(indice),fieldc_user)  
        endrep until match ne -1
        if (match eq 9999) then begin
                print, fieldc_user, ' not found !'
                stop
        endif
;print, 'selected for over-contour:  ',desc_field(indice)
        varcontour=reform(field(indice,*,*,*,*))
endif

;
; overplot wind vectors ?
;
if (n_elements(winds) ge 2) then begin
        ; zonal wind
        match=-1
        indice=-1
        repeat begin
                indice=indice+1
                if (indice eq nfields) $
                        then match=9999 $
                        else match=STRPOS(desc_field(indice),winds(0))  
        endrep until match ne -1
        if (match eq 9999) then begin
                print, 'u not found !'
                stop
        endif
;print, 'selected for vector_x:  ',desc_field(indice)
        u=reform(field(indice,*,*,*,*))
        ; meridional wind
        match=-1
        indice=-1
        repeat begin
                indice=indice+1
                if (indice eq nfields) $
                        then match=9999 $
                        else match=STRPOS(desc_field(indice),winds(1))  
        endrep until match ne -1
        if (match eq 9999) then begin
                print, 'v not found !'
                stop
        endif
;print, 'selected for vector_y:  ',desc_field(indice)
        v=reform(field(indice,*,*,*,*))
endif



;
; loop for graphics   
;
for field_loop=field_loop_start,field_loop_end do begin

        charstart=STRSPLIT(desc_field(field_loop), ' ', LENGTH=charlength)
        charvar=STRMID(desc_field(field_loop),charstart(0),charlength(0))
        var=reform(field(field_loop,*,*,*,*))

        ; useless field ...
        if (charvar eq vert_coord) then goto, skip

        ;
        ; min/max on the whole field
        ;
        w=where(var ne 1e37)
        if (w(0) ne -1) then begin
                minfield=min(var[w])
                maxfield=max(var[w])
                        ;;; manage unvalid points
                        ;;;w=where(var eq 1e37)
                        ;;;if (w(0) ne -1) then var[w]=!Values.F_NAN
        endif else begin
                var=0
                print, 'no valid values in this field', charvar 
        endelse

;----------------------
; initialize graphics
;----------------------

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


lct = round(t[when-1]+24) MOD 24

;;lct = when  ;; on choisit de les numeroter dans l'ordre 
;;if ((when ne -1) and (lct lt 10)) then charend=charend+'0'+string(lct, '(I0)')
;;if (lct ge 10) then charend=charend+string(lct, '(I0)') 
;;lct = round(t[when-1]) MOD 24
;;denom=plot_user+'_'+vert_coord+'_'+charvar
;;if (fieldc_user ne '') then denom=denom+'_'+fieldc_user
;;if (winds(0) ne '') then denom=denom+'_'+winds(0)+winds(1)
;;denom=denom+charend


if (n_elements(file_ps) ne '') then begin
        set_plot, 'ps'
        device, filename=file_ps+'.ps', /color
endif


case charvar of
'HGT': begin
        title='Topography (km)'
        level=1
        pal=16
        format='(F7.0)'
end
'tk': begin
        title='Temperature (K)'
        pal=3
        format='(F4.0)'
end
'TSK': begin
        title='Surface Temperature (K)'
        pal=3
        level=1
end
'U': begin
        title='Zonal wind (!Nm!N.s!U-1!N)'
        pal=33
        if (winds(0) eq 'U') then title='Wind field (!Nm!N.s!U-1!N)'
end
'V': begin
        title='Meridional wind (!Nm!N.s!U-1!N)'
        pal=33
        if (winds(0) eq 'V') then title='Wind field (!Nm!N.s!U-1!N)'
end
'pressure': begin
        title='Pressure (hPa)'
        pal=16
end
'height': begin
        title='Height (km)'
        pal=16
end
'T': title='Potential temperature departure from t0 (K)'
'W': begin
        title='Vertical wind (!Nm!N.s!U-1!N)'
        pal=16
        format='(F6.2)'
end
else: begin
        title=''
        pal=33
end
endcase

if (when ne -1) then begin
        title=title+', LT '+string(lct, '(I0)')
endif else begin
        title=title+', time '+string(when, '(I0)')
endelse
        
title_save=title



;---------------------------------
; choose positioning ...
;---------------------------------

case plot_user of
'map': begin
;
; 1. map: looping on vertical levels
;
        if (level eq -1) then begin
                theloopmin=0
                theloopmax=bottom_top-1
        endif else begin
                theloopmin=level-1 ;!! IDL arrays starts at 0
                theloopmax=level-1
        endelse 
        if (when eq -1) then nt = 0 else nt=when-1
end
'profile': begin
;
; 2. profile: looping on time
;
        ;;discrete=4
        title_axis=strarr(2)
        title_axis(0)=title_save
        if (when eq -1) then begin
                theloopmin=0
                theloopmax=time-1
        endif else begin
                theloopmin=when-1 ;!! IDL arrays starts at 0
                theloopmax=when-1
        endelse 

        column=z
        case vert_coord of
                'height': title_axis(1)='Altitude (km)'
                'pressure': begin
                        H=10.
                        p0=6.1
                        column=H*alog(p0/z)
                        title_axis(1)='Log-Pressure altitude (km)'
                end
                'model_level': title_axis(1)='Model level'
        endcase

        displaylat=float(round(lat(nlat)*100))/100.
        displaylon=float(round(lon(nlon)*100))/100.
        description='WRF Profile at latitude '+string(displaylat,format='(1F6.2)')+$
        '!Uo'+'!N and longitude '+string(displaylon,format='(1F7.2)')+'!Uo!N '
        title=description
end
'hovmuller': begin
;
; 3. hovmuller: looping on vertical levels
;
        if (level eq -1) then begin
                theloopmin=0
                theloopmax=bottom_top-1
        endif else begin
                theloopmin=level-1 ;!! IDL arrays starts at 0
                theloopmax=level-1
        endelse 
end
'zonal': begin
;
; 4. zonal: looping on latitude
;
        if (nlat eq -1) then begin
                theloopmin=0
                theloopmax=south_north-1
        endif else begin
                theloopmin=nlat-1 ;!! IDL arrays starts at 0
                theloopmax=nlat-1
        endelse 
        if (when eq -1) then nt = 0 else nt=when-1
        case vert_coord of
                'height': charlev='Altitude (km)'
                'pressure': begin
                        H=10.
                        p0=6.1
                        z=H*alog(p0/z)
                        charlev='Log-Pressure altitude (km)'
                end
                'model_level': charlev='Model level'
        endcase
end
'meridional': begin
;
; 5. meridional: looping on longitude
;
        if (nlon eq -1) then begin
                theloopmin=0
                theloopmax=west_east-1
        endif else begin
                theloopmin=nlon-1 ;!! IDL arrays starts at 0
                theloopmax=nlon-1
        endelse 
        if (when eq -1) then nt = 0 else nt=when-1
        case vert_coord of
                'height': charlev='Altitude (km)'
                'pressure': begin
                        H=10.
                        p0=6.1
                        z=H*alog(p0/z)
                        charlev='Log-Pressure altitude (km)'
                end
                'model_level': charlev='Model level'
        endcase
end
else: stop
endcase


;-------------
; plot loop
;-------------  
for theloop=theloopmin,theloopmax do begin

case plot_user of
'map': begin
        ;
        ; case 1: map
        ;       
        case vert_coord of
        'height': charlev=string(round(z(theloop)*1000.),'(I0)')+' m'
        'pressure': charlev=string(round(z(theloop)*100.),'(I0)')+' Pa'
        'model_level': charlev=string(round(z(theloop)),'(I0)')
        endcase
        if (n_elements(var(0,0,*,0) ge 2)) then title=title_save+', at level '+charlev
        print, title

        ;-------------
        ; draw map
        ;-------------

        what_I_plot=var[*,*,theloop,nt]

        if (diff_scale eq 0) then begin
                minfield_init=minfield  ;-same scale for all levels 
                maxfield_init=maxfield  ;-same scale for all levels
        endif
        if (fieldc_user eq '') then begin
                overcontour=0.
        endif else begin
                if (topo eq 1) then begin
                        overcontour=varcontour[*,*,0,0]/1000.  ; other levels than 0 are blank !
                endif else begin
                        overcontour=varcontour[*,*,theloop,nt]
                endelse
        endelse
        if (n_elements(winds) lt 2) then begin
                overvector_x=0.
                overvector_y=0.
        endif else begin
                overvector_x=u[*,*,theloop,nt]
                overvector_y=v[*,*,theloop,nt]


;--------------------------------------------------------
; NB: you can use the winds fields to combine variables 
;--------------------------------------------------------
if (n_elements(winds) ge 2) then begin


;** HYDROSTATIC REDUCTION TO LEVEL OF REFERENCE
@sea_level.idl

;** WIND STRESS
@wind_stress.idl

;** ERTEL POTENTIAL VORTICITY
@ertel_vorticity.idl


endif
;--------------------------------------------------------
; NB: you can use the winds fields to combine variables
;--------------------------------------------------------

        endelse

;;******************TEMP
;;pressure anomaly
;nday=3
;oneday=reform(var[*,*,theloop,0:nday-1])
;what_I_plot=what_I_plot-(total(oneday,3)/float(nday))
;;******************TEMP

        map_latlon, $
                what_I_plot, $
                lon, $
                lat, $
                ct=pal, $
                minfield=minfield_init, $
                maxfield=maxfield_init, $
                overcontour=overcontour, $
                overvector_x=overvector_x, $
                overvector_y=overvector_y, $
                colors=colors, $
                title=title, $
                format=format   
end
'profile': begin
        ;
        ; case 2: profile
        ;
        print, title

        what_I_plot=var[nlon,nlat,*,theloop]
        if (diff_scale eq 0) then begin
                minfield_init=minfield  ;-same scale for all levels 
                maxfield_init=maxfield  ;-same scale for all levels
        endif
        if (fieldc_user eq '') then begin
                overcontour=0.
                overcontourz=0.
        endif else begin
                overcontour=varcontour[nlon,nlat,*,theloop]
                overcontourz=0.
        endelse
        if (n_elements(inprofile) gt 1) then begin
                overcontour=inprofile
                overcontourz=0.
                if (n_elements(incolumn) gt 1) then overcontourz=incolumn
        endif

        profile, $
                what_I_plot, $                          ; 1D vertical profile
                column, $                               ; altitudes     
        ;       alt=alt, $                              ; altitude range [altmin, altmax]
                minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
                maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
                inprofile=overcontour, $                ; another vertical profile to overplot
                incolumn=overcontourz, $                ; altitudes of the other vertical profile (in case /= column)
                title_plot=title, $                     ; title of the plot ('Profile' is default)
                title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
        ;       mention=mention, $                      ; add text precision within the plot window (default is nothing or '')
                discrete=discrete                       ; show the profile points (= type of points in !psym)
end
'hovmuller': begin
        ;
        ; case 3: hovmuller
        ;
        ; NB: lat/lon hovmuller 
        ;
        what_I_plot=reform(var[*,nlat,theloop,*])
        space=lon
        temps=findgen(time)*2
        if (fieldc_user ne '') then begin
                overcontour=reform(varcontour[*,nlat,theloop,*])
        endif else begin
                overcontour=0.
        endelse
        minspace=0.
        maxspace=0.

;;******************TEMP
;       minspace=-90.
;       maxspace=90.

if (n_elements(winds) ge 2) then begin          ;;;winds sert à faire des calculs !!
overvector_x=reform(u[*,nlat,4,*])      ;;T 1km
overvector_y=reform(v[*,nlat,0,*])      ;;topo (m)
zref=median(overvector_y)
H=192.*overvector_x/3.72        ;unit is m
;print, (overvector_y-zref)/H
what_I_plot=what_I_plot*exp((overvector_y-zref)/H)


endif 
;oneday=reform(var[*,nlat,0,0:11])
;slice=(total(oneday,2)/12.)
;for i=0, n_elements(var[*,0,0,0])-1 do begin
;       what_I_plot(i,*)=what_I_plot(i,*)-slice(i)
;endfor
;;******************TEMP





        hovmuller, $
                what_I_plot, $                          ; 2D field
                space, $                                ; space coordinate
                temps, $                                ; time coordinate
                minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
                maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
                minspace=minspace, $                    ; minimum value of space window (=0: calculate)
                maxspace=maxspace, $                    ; maximum value of space window (=0: calculate)
                overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=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)
end
'zonal': begin
        ;
        ; case 4: zonal section
        ; 
        displaylat=float(round(lat(theloop)*100))/100.
        title=title_save+ $
                ' section at latitude '+string(displaylat,format='(1F6.2)')+ $
                '!Uo!N '
        print, title

        what_I_plot=reform(var[*,theloop,*,nt])
        space=lon
        altitude=z

if (topo eq 1) then begin
        topography=varcontour[theloop,*,0,0] ; other levels than 0 are blank !
        fieldc_user=''
endif else begin
        topography=0.
endelse 
                    
        
        if (fieldc_user ne '') then begin
                overcontour=reform(varcontour[*,theloop,*,nt])
        endif else begin
                overcontour=0.
        endelse

        if (n_elements(winds) lt 2) then begin
                overvector_x=0.
                overvector_y=0.
        endif else begin
                overvector_x=reform(u[*,theloop,*,nt])          ;;horizontal component
                overvector_y=reform(v[*,theloop,*,nt])          ;;vertical component
        endelse

        minspace=0.
        maxspace=0.

        title_user=title
        title_axis=['Longitude',charlev]

        section, $
                what_I_plot, $                          ; 2D field
                space, $                                ; horizontal coordinate
                altitude, $                             ; altitude coordinate
                minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
                maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
                minspace=minspace, $                    ; minimum value of space window (=0: calculate)
                maxspace=maxspace, $                    ; maximum value of space window (=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)
                colors=colors, $                        ; number of colors/levels (32 is default)
                title_plot=title_user, $                ; title of the plot ('Profile' is default)
                title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
                ct=pal, $                               ; color table (33-rainbow is default)
topo=topography, $
                format=format
end     
'meridional': begin
        ;
        ; case 5: meridional section
        ; 
        displaylon=float(round(lon(theloop)*100))/100.
        title=title_save+ $
                ' section at longitude '+string(displaylon,format='(1F7.2)')+ $
                '!Uo!N '
        print, title

        what_I_plot=reform(var[theloop,*,*,nt])
        space=lat
        altitude=z

if (topo eq 1) then begin
        topography=varcontour[theloop,*,0,0] ; other levels than 0 are blank !
        fieldc_user=''
endif else begin
        topography=0.
endelse
                        
        if (fieldc_user ne '') then begin
                overcontour=reform(varcontour[theloop,*,*,nt])
        endif else begin
                overcontour=0.
        endelse

        if (n_elements(winds) lt 2) then begin
                overvector_x=0.
                overvector_y=0.
        endif else begin
                overvector_x=reform(u[theloop,*,*,nt])          ;;horizontal component
                overvector_y=reform(v[theloop,*,*,nt])          ;;vertical component
        endelse

        minspace=0.
        maxspace=0.

        title_user=title
        title_axis=['Latitude',charlev]

        section, $
                what_I_plot, $                          ; 2D field
                space, $                                ; horizontal coordinate
                altitude, $                             ; altitude coordinate
                minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
                maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
                minspace=minspace, $                    ; minimum value of space window (=0: calculate)
                maxspace=maxspace, $                    ; maximum value of space window (=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)
                colors=colors, $                        ; number of colors/levels (32 is default)
                title_plot=title_user, $                ; title of the plot ('Profile' is default)
                title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
topo=topography, $
                ct=pal, $                               ; color table (33-rainbow is default)
                format=format
end     








else: 
endcase


print, '...ok'
endfor
if (n_elements(file_ps) ne '') then device, /close

        ; save the profiles for another plot (last of the time loop)
        if (n_elements(inprofile) le 1) and (plot_user eq 'profile') then begin
                inprofile=what_I_plot
                incolumn=column
        endif

skip:
endfor

save_data=field
field=0.        ;; free memory

end