Changeset 91 for trunk/mesoscale/PLOT/SPEC/MAP
- Timestamp:
- Mar 8, 2011, 5:59:28 PM (14 years ago)
- Location:
- trunk/mesoscale/PLOT/SPEC/MAP
- Files:
-
- 1 added
- 2 edited
-
defs/map_uvt_inc.pro (modified) (2 diffs)
-
defs/map_uvt_inc.pro.oldref (added)
-
map_uvt.pro (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/mesoscale/PLOT/SPEC/MAP/defs/map_uvt_inc.pro
r85 r91 1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FOLDER 2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FOLDER 3 folder = '/d5/aslmd/LMD_MM_MARS_SIMUS/OM/' 4 coord2d = 'false' ;; for non-regular projections 5 filename = folder + 'OM6_TI85/wrfout_d01_2024-06-43_06:00:00_zabg' 6 save_ps = 'OM6_TI85_tsurf' 7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TIME 8 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TIME 9 ini_utc = 6 ;; cf. name of file 10 freq = 1 ;; cf. 1 output per ... hour 11 utc_to_lt = -9 ;; cf. longitude -- LT = UTC + utc_to_lt 12 use_lt = 2 ;; cf. what user wants 13 use_utc = use_lt - utc_to_lt 14 ntime = (use_utc - ini_utc)/freq ;; TRUE IDL SUBSCRIPT ... ajouter un modulo... 15 print, 'CHECK ... ', filename, ntime 16 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; LVL & FLD 17 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; LVL & FLD 18 nlevel = 0 19 field1 = 'TSURF' ;; comment to trace horizontal velocity 20 no3d = 'true' 21 overvector_x = 0 ;; comment out to get rid of vectors 22 overvector_y = 0 ;; comment out to get rid of vectors 23 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TWEAK VAR 24 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TWEAK VAR 25 ;a0 = 65.052165 & a1 = 3.1228993 & a2 = 0.0053787417 ;; Ls ~ 173 26 ;;a0 = 64.039300 & a1 = 3.1378104 & a2 = 0.0055369148 ;; Ls ~ 120 27 ;what_I_plot = - ( alog ( a1 - what_I_plot / a0 ) ) / a2 28 ;print, max(what_I_plot), min(what_I_plot) 29 ;overvector_x=0 30 ;overvector_y=0 31 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 32 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1 ;;;;;*************************************** FOLDER 2 folder = '/tmp7/aslmd/' 3 coord2d = 'false' ;; for non-regular projections ;;CHANGER windowx and windowy 4 coord2d = 'polar' ;; for polar projection with map_set ;; utiliser use_lt=99 5 filename = folder + 'wrfout_d01_2024-05-02_01:00:00' 6 ini_utc = 1 ;; cf. name of file 7 freq = 1 ;; cf. 1 output per ... hour 8 save_ps = 'POLAR_newphys_RAC_TAUTES' 33 9 10 ;;;;;*************************************** LVL & FLD & TIME 11 nlevel = 0 12 use_lt = 99 ;; cf. what user wants (99 pour tous) 13 field1 = 'TAU_ICE' ;; comment to trace horizontal velocity 14 ;field2 = 'XLAT' ;; contour 15 no3d = 'true' 16 ;overvector_x = 0 ;; comment out to get rid of vectors 17 ;overvector_y = 0 ;; comment out to get rid of vectors 18 19 ;;;;;*************************************** TWEAK VAR 20 ;what_I_plot = what_I_plot * 0. 21 ;print, max(what_I_plot), min(what_I_plot) 34 22 35 23 ;;;;;*************************************** PLOT TITLES 36 ;title_user = 'Surface temperature (K) and Winds 10m ABG (m s!U-1!N)' 37 ;title_user = 'Wind-induced apparent thermal inertia (J m!U-2!N s!U-0.5!N K!U-1!N)' 38 ;title_user = 'Winds 10m ABG (m s!U-1!N)' 39 title_user = 'Surface temperature (K)' 40 ;subtitle_user = 'LT = 03:00am / Ls = 173!Uo!N / dx = 10km [single] / Uniform TI = 85 J m!U-2!N s!U-0.5!N K!U-1!N' 41 subtitle_user = 'LT = 02:00am / Ls = 173!Uo!N / dx = 6km [single] / Uniform TI = 85 J m!U-2!N s!U-0.5!N K!U-1!N' 24 title_user = 'Water ice cloud optical depth at 825 cm!U-1!N' 42 25 title_axis = ['East longitude','North latitude'] 26 ;subtitle_user = 'LMD_MM' 27 ;;subtitle_user = subtitle_user + ' / LT = '+string(use_lt,'(I0)')+'h' 28 ;;subtitle_user = subtitle_user + ' / UTC!D0!N = 06:00am' ;; pour les series 29 ;subtitle_user = subtitle_user + ' / Ls = 120!Uo!N' 30 ;subtitle_user = subtitle_user + ' / dx = 10km [S]' 43 31 44 32 ;;;;;*************************************** COLOR TABLES 45 33 flag_cb = 'true' 46 poscb = 0.80;0.96 47 format = '(I0)' 48 ;format = '(F4.1)' 49 ;colors = 128 34 poscb = 0.85 35 format = '(F5.2)' 50 36 colors = 32 51 37 pal = 22 ;; GOOD: 4, 18, 22, 16, 37, 33, 39, 6, 11 52 38 53 39 ;;;;;*************************************** FILL LIMITS 54 minfield_init = 150. 55 maxfield_init = 180. 56 ;minfield_init = 50. 57 ;maxfield_init = 250. 40 minfield_init = 0. 41 maxfield_init = 0.55 42 ndiv = 11 58 43 59 44 ;;;;;*************************************** LIMIT TRICKS … … 65 50 66 51 ;;;;;*************************************** WINDS 67 ;windex = 10. ;; DEF: 20.68 52 windex = 30. ;; DEF: 20. 69 ;stride = 1. ;; DEF: 5.70 ;stride = 2. ;; DEF: 5.71 53 stride = 3. ;; DEF: 5. 72 54 73 55 ;;;;;*************************************** CONTOUR 74 56 overcontour = overcontour/1000. 75 ;lev = 50. + 50.*findgen(20) 76 ;lev = -10. + 0.2*findgen(20) 77 lev = -10. + 2.*findgen(20) 78 ;lev = -10. + 1.*findgen(40) 57 lev = -10. + 0.5*findgen(80) 79 58 80 59 ;;;;;*************************************** AXIS 81 60 isotropic = 'true' 82 ;intervalx = 0.5 83 intervalx = 2.0 84 ;intervaly = intervalx 85 intervaly = 1.0 61 intervalx = 30.0 62 intervaly = 05.0 86 63 87 64 ;;;;;*************************************** MAP LIMITS 88 ;windowx = [-144.,-126.] 89 ;windowy = [10.,26.] 90 windowx = [-146.,-126.] 91 ;windowy = [10.,28.] 92 ;windowx = [-146.,-128.] 93 ;windowy = [12.,26.] 94 windowy = [11.,27.] 65 windowx = [-180.,180.] 66 windowy = [65.,90.] 67 68 ;;;;;***************************************;;;;; 69 ;;;;;***************************************;;;;; 70 ;;;;;***************************************;;;;; 71 ;;;;;***************************************;;;;; 72 73 ;;;;;*************************************** SETTING TIME (do not modify) 74 if (n_elements(windowx) eq 0) then windowx = 0 75 utc_to_lt = mean(windowx) / 15. ;; cf. longitude -- LT = UTC + utc_to_lt 76 use_utc = use_lt - utc_to_lt 77 zentime = floor(((24 + use_utc - ini_utc) MOD 24)/freq) ;; TRUE IDL SUBSCRIPT... 78 if (use_lt ne 99) then ntime = zentime else ntime = 99 ;; ou commenter pour avoir tous les pas de temps 79 print, zentime, ntime 95 80 96 81 ;;;;;*************************************** -
trunk/mesoscale/PLOT/SPEC/MAP/map_uvt.pro
r86 r91 12 12 ; 13 13 ; 14 zefile=save_ps 14 if (n_elements(field1) ne 0) then getcdf, file=filename, charvar=field1, invar=cfield1 15 if (n_elements(field2) ne 0) then getcdf, file=filename, charvar=field2, invar=cfield2 16 if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin 17 u = getget(filename, 'Um', count=[0,0,1,1], offset=[0,0,nlevel,0]) 18 v = getget(filename, 'Vm', count=[0,0,1,1], offset=[0,0,nlevel,0]) 19 endif 20 getcdf, file=filename, charvar='XLONG', invar=longi 21 getcdf, file=filename, charvar='XLAT', invar=lati 22 getcdf, file=filename, charvar='HGT', invar=hgt 23 ; 24 ; 25 ; 26 sp = 5 ;; relaxation width 27 nx = n_elements(longi(0,*)) 28 ny = n_elements(longi(*,0)) 29 if (n_elements(field1) ne 0) then begin 30 if (no3d ne 'true') then cfield1 = cfield1 (sp:nx-sp-1,sp:ny-sp-1,*,*) else cfield1 = cfield1 (sp:nx-sp-1,sp:ny-sp-1,*) 31 endif 32 if (n_elements(field2) ne 0) then begin 33 if (no3d ne 'true') then cfield2 = cfield2 (sp:nx-sp-1,sp:ny-sp-1,*,*) else cfield2 = cfield2 (sp:nx-sp-1,sp:ny-sp-1,*) 34 endif 35 if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin 36 u = u (sp:nx-sp-1,sp:ny-sp-1,*) 37 v = v (sp:nx-sp-1,sp:ny-sp-1,*) 38 endif 39 longi = longi (sp:nx-sp-1,sp:ny-sp-1,*) 40 lati = lati (sp:nx-sp-1,sp:ny-sp-1,*) 41 hgt = hgt (sp:nx-sp-1,sp:ny-sp-1,*) 42 nx = n_elements(longi(*,0,0)) 43 ny = n_elements(longi(0,*,0)) 44 nt = n_elements(longi(0,0,*)) 45 if ( n_elements( overvector_x ) ne 0 ) then begin 46 overvector_x = u 47 overvector_y = v 48 endif 49 ; 50 ; 51 ; 52 if (ntime eq 99) then begin 53 PRINT, '-- ALL TIME STEPS', nt-1 54 ntstart = 0 & ntend = nt-1 55 endif else begin 56 PRINT, '-- ONLY TIME STEP ', string(ntime,'(I0)') 57 ntstart = ntime & ntend = ntime 58 endelse 59 for ntime = ntstart,ntend do begin 60 ; 61 ; 62 ; 63 zefile=save_ps+string(100+ntime,'(I0)') 15 64 PS_Start, filename=zefile+'.ps' 16 65 print, zefile+'.ps' … … 23 72 ; 24 73 ; 25 if (n_elements(field1) ne 0) then getcdf, file=filename, charvar=field1, invar=cfield126 if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin27 u = getget(filename, 'Um', count=[0,0,1,1], offset=[0,0,nlevel,ntime])28 v = getget(filename, 'Vm', count=[0,0,1,1], offset=[0,0,nlevel,ntime])29 endif30 getcdf, file=filename, charvar='XLONG', invar=longi31 getcdf, file=filename, charvar='XLAT', invar=lati32 getcdf, file=filename, charvar='HGT', invar=hgt33 ;;34 ;;35 ;;36 ; hgt = smooth(hgt, [2,2,0], /EDGE_TRUNCATE) ;;; truc37 ;;38 ;;39 ;;40 sp = 5 ;; relaxation width41 nx = n_elements(longi(0,*))42 ny = n_elements(longi(*,0))43 if (n_elements(field1) ne 0) then begin44 if (no3d ne 'true') then cfield1 = cfield1 (sp:nx-sp-1,sp:ny-sp-1,*,*) else cfield1 = cfield1 (sp:nx-sp-1,sp:ny-sp-1,*)45 endif46 if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin47 u = u (sp:nx-sp-1,sp:ny-sp-1);,*,*)48 v = v (sp:nx-sp-1,sp:ny-sp-1);,*,*)49 endif50 longi = longi (sp:nx-sp-1,sp:ny-sp-1,*)51 lati = lati (sp:nx-sp-1,sp:ny-sp-1,*)52 hgt = hgt (sp:nx-sp-1,sp:ny-sp-1,*)53 nx = n_elements(longi(*,0))54 ny = n_elements(longi(0,*))55 ;56 ;57 ;58 overcontour = reform(hgt(*,*,ntime))59 74 lon = reform(longi(*,*,ntime)) 60 75 lat = reform(lati(*,*,ntime)) 61 if (n_elements(overvector_x) ne 0) then begin62 overvector_x = u;reform(u(*,*,nlevel,ntime))63 overvector_y = v;reform(v(*,*,nlevel,ntime))64 endif65 76 if (n_elements(field1) eq 0) then begin 66 77 print, 'field1: horizontal velocity' … … 69 80 endif else begin 70 81 if (no3d ne 'true') then what_I_plot = reform(cfield1(*,*,nlevel,ntime)) else what_I_plot = reform(cfield1(*,*,ntime)) 82 if (n_elements(u) ne 0) then overvector_x = reform(overvector_x(*,*,ntime)) ;; ne pas utiliser test overvector_x a cause de la boucle temps 83 if (n_elements(v) ne 0) then overvector_y = reform(overvector_y(*,*,ntime)) ;; ne pas utiliser test overvector_y a cause de la boucle temps 84 endelse 85 if (n_elements(field2) eq 0) then begin 86 print, 'field2: topography' 87 overcontour = reform(hgt(*,*,ntime)) 88 endif else begin 89 if (no3d ne 'true') then overcontour = reform(cfield2(*,*,nlevel,ntime)) else overcontour = reform(cfield2(*,*,ntime)) 71 90 endelse 72 91 ; 73 92 ; 74 93 ; 75 minlat=min(lat) & maxlat=max(lat) & minlon=min(lon) & maxlon=max(lon) 76 if (coord2d eq 'true') then begin 94 if ( coord2d eq 'polar' ) then begin 95 print, 'OK YOU USE MAP_SET with POLAR PROJECTION. VECTORS ARE NOT SUPPORTED. USE polar_uv OR ADAPT THIS SCRIPT.' 96 overvector_x = 0 97 overvector_y = 0 98 if (n_elements(windowx) ne 0) then begin 99 latmin = windowy(0) & latmax = windowy(1) & lonmin = windowx(0) & lonmax = windowx(1) 100 endif else begin 101 latmin = 65. & latmax = 90. & lonmin = -180. & lonmax = 180. 102 endelse 103 print, 'latmin,lonmin,latmax,lonmax', latmin,lonmin,latmax,lonmax 104 map_set, latmax, 0., /isotropic, /azimuthal, /noborder, limit=[latmin,lonmin,latmax,lonmax], title=title_user, /advance 105 endif else begin 106 minlat=min(lat) & maxlat=max(lat) & minlon=min(lon) & maxlon=max(lon) 107 if (coord2d eq 'true') then begin 77 108 npoints=n_elements(lon(*,0)) + n_elements(lon(0,*)) ;; trop de points, mais au moins on ne perd rien 78 109 TRIANGULATE, lon, lat, tr 79 110 what_I_plot = GRIDDATA( lon, lat, what_I_plot, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) 80 overvector_x = GRIDDATA( lon, lat, overvector_x, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) 81 overvector_y = GRIDDATA( lon, lat, overvector_y, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) 111 if (n_elements(overvector_x) ne 0) then begin 112 overvector_x = GRIDDATA( lon, lat, overvector_x, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) 113 overvector_y = GRIDDATA( lon, lat, overvector_y, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) 114 endif 82 115 overcontour = GRIDDATA( lon, lat, overcontour, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) 83 ; sale sale sale 84 if (minlat lt min(lat(*,0))) then overvector_y=-overvector_y 85 if (minlon lt min(lon(0,*))) then overvector_x=-overvector_x 116 if (n_elements(overvector_x) ne 0) then begin 117 ; sale sale sale 118 if (minlat lt min(lat(*,0))) then overvector_y=-overvector_y 119 if (minlon lt min(lon(0,*))) then overvector_x=-overvector_x 120 endif 86 121 lon = minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1) 87 122 lat = minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1) 88 endif else begin123 endif else begin 89 124 lon = reform(lon(*,0)) 90 125 lat = reform(lat(0,*)) … … 96 131 ;lon = minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1) 97 132 ;lat = minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1) 133 endelse 98 134 endelse 99 135 help, what_I_plot, lon, lat 100 136 print, min(what_I_plot) 101 137 print, max(what_I_plot) 138 print, min(overcontour) 139 print, max(overcontour) 102 140 ; 103 141 ; … … 119 157 ; 120 158 ; 159 if ( coord2d eq 'polar' ) then begin ;;; pourrait aller dans map_latlon 160 loadct, 0 161 MAP_GRID, $ 162 CHARSIZE = 1., $ 163 COLOR = 0, $ 164 ;lats=-60, $ 165 LABEL = 1, $ ;; /LABEL or LABEL=2 (one label any 2 grid lines) 166 LATDEL = intervaly, $ ;;5 10 167 LONDEL = intervalx, $ ;;15 168 LONLAB = latmin + intervaly/2., $ ;5. + (latmin+latmax)/2., $ ;0. 169 LATLAB = -0.001, $ 170 GLINESTYLE = 2, $ 171 GLINETHICK = 0.3 172 ;LONALIGN = 0., $ 173 ;LATALIGN = 1. 174 endif 175 ; 176 ; 177 ; 121 178 PS_End, /PNG 179 endfor 122 180 end
Note: See TracChangeset
for help on using the changeset viewer.
