Index: /trunk/000-USERS =================================================================== --- /trunk/000-USERS (revision 90) +++ /trunk/000-USERS (revision 91) @@ -19,2 +19,3 @@ Francois Forget LMD [CNRS] 03 - 03 - 2011 All GCMs Jean-Baptiste Madeleine LMD [CNRS] 03 - 03 - 2011 Mars GCM +Julien Faure LMD [CNRS/UPMC] 07 - 03 - 2011 Mars GCM, MMM and LES Index: /trunk/mars/libf/phymars/initracer.F =================================================================== --- /trunk/mars/libf/phymars/initracer.F (revision 90) +++ /trunk/mars/libf/phymars/initracer.F (revision 91) @@ -46,8 +46,9 @@ #include "chimiedata.h" +#include "surfdat.h" real qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2) real co2ice(ngridmx) ! co2 ice mass on surface (e.g. kg.m-2) - integer iq,ig,count + integer iq,ig,count, yeyey real r0_lift , reff_lift, nueff_lift c Ratio of small over large dust particles (used when both @@ -501,5 +502,42 @@ c Perennial H20 north cap defined by watercaptag=true (allows surface to be c hollowed by sublimation in vdifc). + yeyey = 0 do ig=1,ngridmx +! !!! TESTS TESTS outliers +! !!! TESTS TESTS outliers +! if ( ( lati(ig)*180./pi .ge. 75 ) .and. +! . ( lati(ig)*180./pi .le. 77 ) .and. +! . ( ( ( long(ig)*180./pi .ge. 000. ) .and. +! . ( long(ig)*180./pi .le. 120. ) ) +! . .or. +! . ( ( long(ig)*180./pi .ge. -130. ) .and. +! . ( long(ig)*180./pi .le. -115. ) ) ) ) then +! if (yeyey .eq. 0) then !!! 1/2 en 64x48 sinon trop large en lat +! write(*,*) "outliers ", lati(ig)*180./pi, long(ig)*180./pi +! if (ngridmx.ne.1) watercaptag(ig)=.true. +! dryness(ig) = 1. +! albedodat(ig) = 0.45 !! comme alb_surfice +! yeyey = 1 +! else +! yeyey = 0 +! endif +! endif +! !!! TESTS TESTS outliers +! !!! TESTS TESTS outliers +! +! !!! TESTS TESTS addcap +! !!! TESTS TESTS addcap +! if ( ( lati(ig)*180./pi .ge. 82 ) .and. +! . ( lati(ig)*180./pi .le. 84 ) .and. +! . ( ( long(ig)*180./pi .gt. -030. ) .and. +! . ( long(ig)*180./pi .lt. 090. ) ) ) then +! write(*,*) "capadd ", lati(ig)*180./pi, long(ig)*180./pi +! if (ngridmx.ne.1) watercaptag(ig)=.true. +! albedodat(ig) = 0.45 !! comme alb_surfice +! dryness(ig) = 1. +! endif +! !!! TESTS TESTS addcap +! !!! TESTS TESTS addcap + if (lati(ig)*180./pi.gt.84) then if (ngridmx.ne.1) watercaptag(ig)=.true. Index: /trunk/mesoscale/LMDZ.MARS.new/libf/phymars/physiq.F =================================================================== --- /trunk/mesoscale/LMDZ.MARS.new/libf/phymars/physiq.F (revision 90) +++ /trunk/mesoscale/LMDZ.MARS.new/libf/phymars/physiq.F (revision 91) @@ -1338,10 +1338,10 @@ & qsurf(1:ngridmx,igcm_h2o_ice)) -c CALL WRITEDIAGFI(ngridmx,'mtot', -c & 'total mass of water vapor', -c & 'kg/m2',2,mtot) -c CALL WRITEDIAGFI(ngridmx,'icetot', -c & 'total mass of water ice', -c & 'kg/m2',2,icetot) + CALL WRITEDIAGFI(ngridmx,'mtot', + & 'total mass of water vapor', + & 'kg/m2',2,mtot) + CALL WRITEDIAGFI(ngridmx,'icetot', + & 'total mass of water ice', + & 'kg/m2',2,icetot) cc vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) cc & *mugaz/mmol(igcm_h2o_ice) Index: /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F =================================================================== --- /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F (revision 90) +++ /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F (revision 91) @@ -334,5 +334,4 @@ logical tdiag - REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx) real hco2(nqmx),tmean, zlocal(nlayermx) @@ -341,4 +340,9 @@ REAL time_phys + +!!! WRF for retrocompatibility with newphys + REAL tauTES(ngridmx) ! column optical depth at 825 cm-1 + REAL qsurfice(ngridmx) ! pour diagnostics +!!! WRF for retrocompatibility with newphys c======================================================================= Index: /trunk/mesoscale/NOTES.txt =================================================================== --- /trunk/mesoscale/NOTES.txt (revision 90) +++ /trunk/mesoscale/NOTES.txt (revision 91) @@ -1,10 +1,9 @@ +concat.e puis localtime.e puis +ncwa -O -v mtot,icetot -a longitude -d longitude,-179.0,179.0 diagfi.nc yeye.nc +ncwa -O -v mtot -a longitude -d longitude,-180.0,180.0 concat_LT.nc mawd.nc -ncwa -O -v mtot,icetot -a longitude -d longitude,-179.0,179.0 diagfi.nc yeye.nc - - - -A FAIRE:::: mettre des flags precompilo dans les meso_ -les reporter dans makegcm + A FAIRE:::: mettre des flags precompilo dans les meso_ + les reporter dans makegcm changer le renormalisation dans aeropacity ???? Index: /trunk/mesoscale/PLOT/MINIMAL/map_latlon.pro =================================================================== --- /trunk/mesoscale/PLOT/MINIMAL/map_latlon.pro (revision 90) +++ /trunk/mesoscale/PLOT/MINIMAL/map_latlon.pro (revision 91) @@ -73,4 +73,5 @@ endelse if (n_elements(isotropic) eq 0) then isotropic='true' +if (n_elements(coord2d) eq 0) then coord2d='false' ;------------------ @@ -160,5 +161,5 @@ if (isotropic eq 'true') then begin loadct, 0, /silent - if ( n_elements(polar) eq 0 ) then begin + if ( coord2d ne 'polar' ) then begin ;;trace non polaire contour, /nodata, $ @@ -183,9 +184,8 @@ endelse subtitle_user = '' - print, 'YORGL YORGL' endif else begin pospos = [0.15, 0.35, 0.95, 0.95] ;; non-iso loadct, 0, /silent - if ( n_elements(polar) eq 0 ) then begin + if ( coord2d ne 'polar' ) then begin ;;trace non polaire contour, /nodata, $ Index: /trunk/mesoscale/PLOT/SPEC/LES/turb_section_period.pro =================================================================== --- /trunk/mesoscale/PLOT/SPEC/LES/turb_section_period.pro (revision 91) +++ /trunk/mesoscale/PLOT/SPEC/LES/turb_section_period.pro (revision 91) @@ -0,0 +1,53 @@ +pro turb_section_period +; +; +; +what_I_plot=0. +overcontour=0. +SPAWN, '\rm param_plot.idl ; cp turb_section_period_inc.pro param_plot.idl' +; +; +; +filefile = "./MERIDIANI_tau050_wind10_section_plume.nc" +getcdf, file=filefile, charvar="W", invar=zew +getcdf, file=filefile, charvar="PHTOT", invar=zevar +zew = reform(zew) +zevar = reform(zevar) + + nx = n_elements(zevar(*,0)) + nz = n_elements(zevar(0,*)) + + xx = findgen(nx) * 100. / 1000. + zz = reform(zevar(0,*)) / 3.72 / 1000. + ( 1473.35 / 1000. ) + + zew = zew[*,0:nz-2] + zz = zz[0:nz-2] + + + PS_Start, filename='sectionLES.ps' + !P.Charsize = 1.2 + !p.charthick = 2.0 + !p.thick = 2.0 + !x.thick = 2.0 + !y.thick = 2.0 + + what_I_plot = zew + lon = xx + lat = zz + + map_latlon, $ + what_I_plot, $ ; 2D field + lon, $ ; 1D latitude OR 2D + lat, $ ; 1D longitude OR 2D + 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_user, $ ; title of the plot ('' is default) + format=format ; format of colorbar annotations ('(F6.2)' is default) + +PS_End, /PNG +end Index: /trunk/mesoscale/PLOT/SPEC/LES/turb_section_period_inc.pro =================================================================== --- /trunk/mesoscale/PLOT/SPEC/LES/turb_section_period_inc.pro (revision 91) +++ /trunk/mesoscale/PLOT/SPEC/LES/turb_section_period_inc.pro (revision 91) @@ -0,0 +1,22 @@ +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +format = '(I0)' +colors = 64 ;32 +minfield_init = -8. +maxfield_init = 12. +pal = 33 +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +missing_value=1.e30 +lim_max = maxfield_init & w=where((what_I_plot ge lim_max) and (what_I_plot le 1e9)) & if (w[0] ne -1) then what_I_plot[w]=lim_max +lim_min = minfield_init & w=where(what_I_plot le lim_min) & if (w[0] ne -1) then what_I_plot[w]=lim_min +flag_cb = 'true' +windowx = [0.,9.] +windowy = [0.,12.] +intervalx = 1. +intervaly = 1. +title_axis = ['Horizontal coordinate (km)','Altitude above surface (km)'] +title_user = 'Vertical wind (m s!U-1!N) in Meridiani LES' +subtitle_user = 'DOD = 0.5 / BW = 10 m s!U-1!N' +poscb = 0.55 +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Index: /trunk/mesoscale/PLOT/SPEC/MAP/defs/map_uvt_inc.pro =================================================================== --- /trunk/mesoscale/PLOT/SPEC/MAP/defs/map_uvt_inc.pro (revision 90) +++ /trunk/mesoscale/PLOT/SPEC/MAP/defs/map_uvt_inc.pro (revision 91) @@ -1,59 +1,44 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FOLDER -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FOLDER - folder = '/d5/aslmd/LMD_MM_MARS_SIMUS/OM/' - coord2d = 'false' ;; for non-regular projections -filename = folder + 'OM6_TI85/wrfout_d01_2024-06-43_06:00:00_zabg' -save_ps = 'OM6_TI85_tsurf' -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TIME -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TIME -ini_utc = 6 ;; cf. name of file -freq = 1 ;; cf. 1 output per ... hour -utc_to_lt = -9 ;; cf. longitude -- LT = UTC + utc_to_lt - use_lt = 2 ;; cf. what user wants - use_utc = use_lt - utc_to_lt - ntime = (use_utc - ini_utc)/freq ;; TRUE IDL SUBSCRIPT ... ajouter un modulo... - print, 'CHECK ... ', filename, ntime -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; LVL & FLD -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; LVL & FLD -nlevel = 0 -field1 = 'TSURF' ;; comment to trace horizontal velocity -no3d = 'true' -overvector_x = 0 ;; comment out to get rid of vectors -overvector_y = 0 ;; comment out to get rid of vectors -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TWEAK VAR -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TWEAK VAR - ;a0 = 65.052165 & a1 = 3.1228993 & a2 = 0.0053787417 ;; Ls ~ 173 - ;;a0 = 64.039300 & a1 = 3.1378104 & a2 = 0.0055369148 ;; Ls ~ 120 - ;what_I_plot = - ( alog ( a1 - what_I_plot / a0 ) ) / a2 - ;print, max(what_I_plot), min(what_I_plot) - ;overvector_x=0 - ;overvector_y=0 -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;*************************************** FOLDER + folder = '/tmp7/aslmd/' + coord2d = 'false' ;; for non-regular projections ;;CHANGER windowx and windowy + coord2d = 'polar' ;; for polar projection with map_set ;; utiliser use_lt=99 + filename = folder + 'wrfout_d01_2024-05-02_01:00:00' + ini_utc = 1 ;; cf. name of file + freq = 1 ;; cf. 1 output per ... hour + save_ps = 'POLAR_newphys_RAC_TAUTES' +;;;;;*************************************** LVL & FLD & TIME + nlevel = 0 + use_lt = 99 ;; cf. what user wants (99 pour tous) + field1 = 'TAU_ICE' ;; comment to trace horizontal velocity + ;field2 = 'XLAT' ;; contour + no3d = 'true' + ;overvector_x = 0 ;; comment out to get rid of vectors + ;overvector_y = 0 ;; comment out to get rid of vectors + +;;;;;*************************************** TWEAK VAR + ;what_I_plot = what_I_plot * 0. + ;print, max(what_I_plot), min(what_I_plot) ;;;;;*************************************** PLOT TITLES - ;title_user = 'Surface temperature (K) and Winds 10m ABG (m s!U-1!N)' - ;title_user = 'Wind-induced apparent thermal inertia (J m!U-2!N s!U-0.5!N K!U-1!N)' - ;title_user = 'Winds 10m ABG (m s!U-1!N)' - title_user = 'Surface temperature (K)' - ;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' - 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' + title_user = 'Water ice cloud optical depth at 825 cm!U-1!N' title_axis = ['East longitude','North latitude'] + ;subtitle_user = 'LMD_MM' + ;;subtitle_user = subtitle_user + ' / LT = '+string(use_lt,'(I0)')+'h' + ;;subtitle_user = subtitle_user + ' / UTC!D0!N = 06:00am' ;; pour les series + ;subtitle_user = subtitle_user + ' / Ls = 120!Uo!N' + ;subtitle_user = subtitle_user + ' / dx = 10km [S]' ;;;;;*************************************** COLOR TABLES flag_cb = 'true' - poscb = 0.80;0.96 - format = '(I0)' - ;format = '(F4.1)' - ;colors = 128 + poscb = 0.85 + format = '(F5.2)' colors = 32 pal = 22 ;; GOOD: 4, 18, 22, 16, 37, 33, 39, 6, 11 ;;;;;*************************************** FILL LIMITS - minfield_init = 150. - maxfield_init = 180. - ;minfield_init = 50. - ;maxfield_init = 250. + minfield_init = 0. + maxfield_init = 0.55 + ndiv = 11 ;;;;;*************************************** LIMIT TRICKS @@ -65,32 +50,32 @@ ;;;;;*************************************** WINDS - ;windex = 10. ;; DEF: 20. windex = 30. ;; DEF: 20. - ;stride = 1. ;; DEF: 5. - ;stride = 2. ;; DEF: 5. stride = 3. ;; DEF: 5. ;;;;;*************************************** CONTOUR overcontour = overcontour/1000. - ;lev = 50. + 50.*findgen(20) - ;lev = -10. + 0.2*findgen(20) - lev = -10. + 2.*findgen(20) - ;lev = -10. + 1.*findgen(40) + lev = -10. + 0.5*findgen(80) ;;;;;*************************************** AXIS isotropic = 'true' - ;intervalx = 0.5 - intervalx = 2.0 - ;intervaly = intervalx - intervaly = 1.0 + intervalx = 30.0 + intervaly = 05.0 ;;;;;*************************************** MAP LIMITS - ;windowx = [-144.,-126.] - ;windowy = [10.,26.] - windowx = [-146.,-126.] - ;windowy = [10.,28.] - ;windowx = [-146.,-128.] - ;windowy = [12.,26.] - windowy = [11.,27.] + windowx = [-180.,180.] + windowy = [65.,90.] + +;;;;;***************************************;;;;; +;;;;;***************************************;;;;; +;;;;;***************************************;;;;; +;;;;;***************************************;;;;; + +;;;;;*************************************** SETTING TIME (do not modify) + if (n_elements(windowx) eq 0) then windowx = 0 + utc_to_lt = mean(windowx) / 15. ;; cf. longitude -- LT = UTC + utc_to_lt + use_utc = use_lt - utc_to_lt + zentime = floor(((24 + use_utc - ini_utc) MOD 24)/freq) ;; TRUE IDL SUBSCRIPT... + if (use_lt ne 99) then ntime = zentime else ntime = 99 ;; ou commenter pour avoir tous les pas de temps + print, zentime, ntime ;;;;;*************************************** Index: /trunk/mesoscale/PLOT/SPEC/MAP/defs/map_uvt_inc.pro.oldref =================================================================== --- /trunk/mesoscale/PLOT/SPEC/MAP/defs/map_uvt_inc.pro.oldref (revision 91) +++ /trunk/mesoscale/PLOT/SPEC/MAP/defs/map_uvt_inc.pro.oldref (revision 91) @@ -0,0 +1,139 @@ +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FOLDER +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FOLDER + folder = '/d5/aslmd/LMD_MM_MARS_SIMUS/OM/' + coord2d = 'false' ;; for non-regular projections +filename = folder + 'OM6_TI85/wrfout_d01_2024-06-43_06:00:00_zabg' +save_ps = 'OM6_TI85_tsurf' +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TIME +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TIME +ini_utc = 6 ;; cf. name of file +freq = 1 ;; cf. 1 output per ... hour +utc_to_lt = -9 ;; cf. longitude -- LT = UTC + utc_to_lt + use_lt = 2 ;; cf. what user wants + use_utc = use_lt - utc_to_lt + ntime = (use_utc - ini_utc)/freq ;; TRUE IDL SUBSCRIPT ... ajouter un modulo... + print, 'CHECK ... ', filename, ntime +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; LVL & FLD +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; LVL & FLD +nlevel = 0 +field1 = 'TSURF' ;; comment to trace horizontal velocity +no3d = 'true' +overvector_x = 0 ;; comment out to get rid of vectors +overvector_y = 0 ;; comment out to get rid of vectors +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TWEAK VAR +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; TWEAK VAR + ;a0 = 65.052165 & a1 = 3.1228993 & a2 = 0.0053787417 ;; Ls ~ 173 + ;;a0 = 64.039300 & a1 = 3.1378104 & a2 = 0.0055369148 ;; Ls ~ 120 + ;what_I_plot = - ( alog ( a1 - what_I_plot / a0 ) ) / a2 + ;print, max(what_I_plot), min(what_I_plot) + ;overvector_x=0 + ;overvector_y=0 +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + + +;;;;;*************************************** PLOT TITLES + ;title_user = 'Surface temperature (K) and Winds 10m ABG (m s!U-1!N)' + ;title_user = 'Wind-induced apparent thermal inertia (J m!U-2!N s!U-0.5!N K!U-1!N)' + ;title_user = 'Winds 10m ABG (m s!U-1!N)' + title_user = 'Surface temperature (K)' + ;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' + 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' + title_axis = ['East longitude','North latitude'] + +;;;;;*************************************** COLOR TABLES + flag_cb = 'true' + poscb = 0.80;0.96 + format = '(I0)' + ;format = '(F4.1)' + ;colors = 128 + colors = 32 + pal = 22 ;; GOOD: 4, 18, 22, 16, 37, 33, 39, 6, 11 + +;;;;;*************************************** FILL LIMITS + minfield_init = 150. + maxfield_init = 180. + ;minfield_init = 50. + ;maxfield_init = 250. + + ;;;;;*************************************** LIMIT TRICKS + ;;;;;********************* must always follow FILL LIMITS + missing_value=1.e30 + lim_max = maxfield_init & w=where((what_I_plot ge lim_max) and (what_I_plot le 1e9)) & if (w[0] ne -1) then what_I_plot[w]=lim_max + lim_min = minfield_init & w=where(what_I_plot le lim_min) & if (w[0] ne -1) then what_I_plot[w]=lim_min + ;lim_blank = 2. & w=where(abs(what_I_plot) le lim_blank) & if (w[0] ne -1) then what_I_plot[w]=missing_value + +;;;;;*************************************** WINDS + ;windex = 10. ;; DEF: 20. + windex = 30. ;; DEF: 20. + ;stride = 1. ;; DEF: 5. + ;stride = 2. ;; DEF: 5. + stride = 3. ;; DEF: 5. + +;;;;;*************************************** CONTOUR + overcontour = overcontour/1000. + ;lev = 50. + 50.*findgen(20) + ;lev = -10. + 0.2*findgen(20) + lev = -10. + 2.*findgen(20) + ;lev = -10. + 1.*findgen(40) + +;;;;;*************************************** AXIS + isotropic = 'true' + ;intervalx = 0.5 + intervalx = 2.0 + ;intervaly = intervalx + intervaly = 1.0 + +;;;;;*************************************** MAP LIMITS + ;windowx = [-144.,-126.] + ;windowy = [10.,26.] + windowx = [-146.,-126.] + ;windowy = [10.,28.] + ;windowx = [-146.,-128.] + ;windowy = [12.,26.] + windowy = [11.,27.] + +;;;;;*************************************** +;;;;;*************************************** +;;;;;*************************************** +;;;;;*************************************** + +;;; +;;; VERTICAL SECTION LIMITS +;;; +;minalt=-5. ;; grepSEC +;maxalt=40. ;; grepSEC +;minspace=0.0 ;; grepSEC +;maxspace=35.0 ;; grepSEC + +;;; +;;; METRIC UNITS FOR VERTICAL SECTION +;;; +;factor=10. ;; grepSEC +;space=space*60. & spacekm='true' ;; grepSEC +;minspace=minspace*60./factor ;; grepSEC +;maxspace=maxspace*60./factor ;; grepSEC +;intervalx=round(intervalx*60./factor) ;; grepSEC + + +;;; +;;; TRICKS +;;; + +;; pour tracer juste les vecteurs sur un fond uni ou vide +;pal=0 ;; 1/4 grepMAP +;what_I_plot(*,*)=what_I_plot(*,*)*0.+0.2 ;; 2/4 grepMAP +;what_I_plot(0,0)=0. ;; 3/4 grepMAP +;flag_cb='false' ;; 4/4 grepMAP + +;;; truc pour tracer juste les contours et une zone grisee de topo +;pal=0 ;; 1/9 grepALL +;w=where(abs(what_I_plot) lt missing_value) ;; 2/9 grepALL +;what_I_plot[w]=0. ;; 3/9 grepALL +;w=where(abs(what_I_plot) gt missing_value) ;; 4/9 grepALL +;what_I_plot[w]=0.75 ;; 5/9 grepALL +;w=where(what_I_plot eq 0.) ;; 6/9 grepALL +;what_I_plot[w]=missing_value ;; 7/9 grepALL +;what_I_plot(0,0)=1. & what_I_plot(1,0)=0. ;; 8/9 grepALL +;flag_cb='false' ;; 9/9 grepALL + Index: /trunk/mesoscale/PLOT/SPEC/MAP/map_uvt.pro =================================================================== --- /trunk/mesoscale/PLOT/SPEC/MAP/map_uvt.pro (revision 90) +++ /trunk/mesoscale/PLOT/SPEC/MAP/map_uvt.pro (revision 91) @@ -12,5 +12,54 @@ ; ; - zefile=save_ps +if (n_elements(field1) ne 0) then getcdf, file=filename, charvar=field1, invar=cfield1 +if (n_elements(field2) ne 0) then getcdf, file=filename, charvar=field2, invar=cfield2 +if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin + u = getget(filename, 'Um', count=[0,0,1,1], offset=[0,0,nlevel,0]) + v = getget(filename, 'Vm', count=[0,0,1,1], offset=[0,0,nlevel,0]) +endif +getcdf, file=filename, charvar='XLONG', invar=longi +getcdf, file=filename, charvar='XLAT', invar=lati +getcdf, file=filename, charvar='HGT', invar=hgt +; +; +; +sp = 5 ;; relaxation width +nx = n_elements(longi(0,*)) +ny = n_elements(longi(*,0)) +if (n_elements(field1) ne 0) then begin + 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,*) +endif +if (n_elements(field2) ne 0) then begin + 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,*) +endif +if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin + u = u (sp:nx-sp-1,sp:ny-sp-1,*) + v = v (sp:nx-sp-1,sp:ny-sp-1,*) +endif +longi = longi (sp:nx-sp-1,sp:ny-sp-1,*) +lati = lati (sp:nx-sp-1,sp:ny-sp-1,*) +hgt = hgt (sp:nx-sp-1,sp:ny-sp-1,*) +nx = n_elements(longi(*,0,0)) +ny = n_elements(longi(0,*,0)) +nt = n_elements(longi(0,0,*)) +if ( n_elements( overvector_x ) ne 0 ) then begin + overvector_x = u + overvector_y = v +endif +; +; +; + if (ntime eq 99) then begin + PRINT, '-- ALL TIME STEPS', nt-1 + ntstart = 0 & ntend = nt-1 + endif else begin + PRINT, '-- ONLY TIME STEP ', string(ntime,'(I0)') + ntstart = ntime & ntend = ntime + endelse + for ntime = ntstart,ntend do begin +; +; +; + zefile=save_ps+string(100+ntime,'(I0)') PS_Start, filename=zefile+'.ps' print, zefile+'.ps' @@ -23,44 +72,6 @@ ; ; -if (n_elements(field1) ne 0) then getcdf, file=filename, charvar=field1, invar=cfield1 -if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin - u = getget(filename, 'Um', count=[0,0,1,1], offset=[0,0,nlevel,ntime]) - v = getget(filename, 'Vm', count=[0,0,1,1], offset=[0,0,nlevel,ntime]) -endif -getcdf, file=filename, charvar='XLONG', invar=longi -getcdf, file=filename, charvar='XLAT', invar=lati -getcdf, file=filename, charvar='HGT', invar=hgt -;; -;; -;; -; hgt = smooth(hgt, [2,2,0], /EDGE_TRUNCATE) ;;; truc -;; -;; -;; -sp = 5 ;; relaxation width -nx = n_elements(longi(0,*)) -ny = n_elements(longi(*,0)) -if (n_elements(field1) ne 0) then begin - 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,*) -endif -if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin - u = u (sp:nx-sp-1,sp:ny-sp-1);,*,*) - v = v (sp:nx-sp-1,sp:ny-sp-1);,*,*) -endif -longi = longi (sp:nx-sp-1,sp:ny-sp-1,*) -lati = lati (sp:nx-sp-1,sp:ny-sp-1,*) -hgt = hgt (sp:nx-sp-1,sp:ny-sp-1,*) -nx = n_elements(longi(*,0)) -ny = n_elements(longi(0,*)) -; -; -; -overcontour = reform(hgt(*,*,ntime)) lon = reform(longi(*,*,ntime)) lat = reform(lati(*,*,ntime)) -if (n_elements(overvector_x) ne 0) then begin - overvector_x = u;reform(u(*,*,nlevel,ntime)) - overvector_y = v;reform(v(*,*,nlevel,ntime)) -endif if (n_elements(field1) eq 0) then begin print, 'field1: horizontal velocity' @@ -69,22 +80,46 @@ endif else begin if (no3d ne 'true') then what_I_plot = reform(cfield1(*,*,nlevel,ntime)) else what_I_plot = reform(cfield1(*,*,ntime)) + 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 + 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 +endelse +if (n_elements(field2) eq 0) then begin + print, 'field2: topography' + overcontour = reform(hgt(*,*,ntime)) +endif else begin + if (no3d ne 'true') then overcontour = reform(cfield2(*,*,nlevel,ntime)) else overcontour = reform(cfield2(*,*,ntime)) endelse ; ; ; -minlat=min(lat) & maxlat=max(lat) & minlon=min(lon) & maxlon=max(lon) -if (coord2d eq 'true') then begin +if ( coord2d eq 'polar' ) then begin + print, 'OK YOU USE MAP_SET with POLAR PROJECTION. VECTORS ARE NOT SUPPORTED. USE polar_uv OR ADAPT THIS SCRIPT.' + overvector_x = 0 + overvector_y = 0 + if (n_elements(windowx) ne 0) then begin + latmin = windowy(0) & latmax = windowy(1) & lonmin = windowx(0) & lonmax = windowx(1) + endif else begin + latmin = 65. & latmax = 90. & lonmin = -180. & lonmax = 180. + endelse + print, 'latmin,lonmin,latmax,lonmax', latmin,lonmin,latmax,lonmax + map_set, latmax, 0., /isotropic, /azimuthal, /noborder, limit=[latmin,lonmin,latmax,lonmax], title=title_user, /advance +endif else begin + minlat=min(lat) & maxlat=max(lat) & minlon=min(lon) & maxlon=max(lon) + if (coord2d eq 'true') then begin npoints=n_elements(lon(*,0)) + n_elements(lon(0,*)) ;; trop de points, mais au moins on ne perd rien TRIANGULATE, lon, lat, tr what_I_plot = GRIDDATA( lon, lat, what_I_plot, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) - overvector_x = GRIDDATA( lon, lat, overvector_x, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) - overvector_y = GRIDDATA( lon, lat, overvector_y, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) + if (n_elements(overvector_x) ne 0) then begin + overvector_x = GRIDDATA( lon, lat, overvector_x, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) + overvector_y = GRIDDATA( lon, lat, overvector_y, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) + endif overcontour = GRIDDATA( lon, lat, overcontour, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN ) - ; sale sale sale - if (minlat lt min(lat(*,0))) then overvector_y=-overvector_y - if (minlon lt min(lon(0,*))) then overvector_x=-overvector_x + if (n_elements(overvector_x) ne 0) then begin + ; sale sale sale + if (minlat lt min(lat(*,0))) then overvector_y=-overvector_y + if (minlon lt min(lon(0,*))) then overvector_x=-overvector_x + endif lon = minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1) lat = minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1) -endif else begin + endif else begin lon = reform(lon(*,0)) lat = reform(lat(0,*)) @@ -96,8 +131,11 @@ ;lon = minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1) ;lat = minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1) + endelse endelse help, what_I_plot, lon, lat print, min(what_I_plot) print, max(what_I_plot) +print, min(overcontour) +print, max(overcontour) ; ; @@ -119,4 +157,24 @@ ; ; +if ( coord2d eq 'polar' ) then begin ;;; pourrait aller dans map_latlon + loadct, 0 + MAP_GRID, $ + CHARSIZE = 1., $ + COLOR = 0, $ + ;lats=-60, $ + LABEL = 1, $ ;; /LABEL or LABEL=2 (one label any 2 grid lines) + LATDEL = intervaly, $ ;;5 10 + LONDEL = intervalx, $ ;;15 + LONLAB = latmin + intervaly/2., $ ;5. + (latmin+latmax)/2., $ ;0. + LATLAB = -0.001, $ + GLINESTYLE = 2, $ + GLINETHICK = 0.3 + ;LONALIGN = 0., $ + ;LATALIGN = 1. +endif +; +; +; PS_End, /PNG + endfor end Index: /trunk/mesoscale/PLOT/SPEC/POLAR/polar.README =================================================================== --- /trunk/mesoscale/PLOT/SPEC/POLAR/polar.README (revision 90) +++ /trunk/mesoscale/PLOT/SPEC/POLAR/polar.README (revision 91) @@ -5,5 +5,5 @@ profile_uvt pour un profil -map_uvt est un vieux truc +map_uvt est un vieux truc, non OK pour tracer hors des poles Index: /trunk/mesoscale/PLOT/SPEC/SERIES/mawd.pro =================================================================== --- /trunk/mesoscale/PLOT/SPEC/SERIES/mawd.pro (revision 91) +++ /trunk/mesoscale/PLOT/SPEC/SERIES/mawd.pro (revision 91) @@ -0,0 +1,91 @@ +pro mawd + +;;;;;;; +f_user='/donnees/aslmd/SVN/trunk/mesoscale/TMPDIR/GCMINI/mawd.nc' +;f_user='/tmp7/aslmd/mawd.nc' +v_user = 'mtot' +x_user = 'Time' +y_user = 'latitude' +transp = 'yes' +tes = 'no' +;;;;;;; +f_user='/donnees/aslmd/TES/TES.SeasonalClimatology.nc' +v_user = 'water' +x_user = 'time' +y_user = 'latitude' +transp = 'yes' +tes = 'yes' +;;;;;;; + + +what_I_plot=0. +overcontour=0 +@mawd_inc.pro +print, lev +SPAWN, '\rm param_plot.idl ; cp mawd_inc.pro param_plot.idl' + +PS_START, file='mawd.ps' + !P.Charsize = 1.2 + !p.charthick = 2.0 + !p.thick = 2.0 + !x.thick = 2.0 + !y.thick = 2.0 + +cdfid = ncdf_open(f_user) + +varid=ncdf_varid(cdfid,v_user) +ncdf_varget, cdfid, varid, champ + +varid=ncdf_varid(cdfid,x_user) +ncdf_varget, cdfid, varid, xx + +varid=ncdf_varid(cdfid,y_user) +ncdf_varget, cdfid, varid, yy + +if (tes eq 'no') then begin + ;;; en precip-microns + champ = champ * 1.e6 / 917. + ;;; entre 0 et 360 + xx = xx MOD 360 +endif else begin + xx = xx - 360.*2. +endelse + +help, champ +if (tes eq 'no') then begin + what_I_plot = champ +endif else begin + what_I_plot = float(reform(champ(0,*,*))) +endelse +help, what_I_plot +overcontour = what_I_plot + +if (tes eq 'yes') then begin + w=where(what_I_plot eq -1.) + what_I_plot[w] = !VALUES.F_NAN + overcontour[w] = !VALUES.F_NAN +endif + +if (transp eq 'yes') then begin + what_I_plot = transpose(what_I_plot) + overcontour = transpose(overcontour) +endif + +map_latlon, $ + what_I_plot, $ ; 2D field + xx, $ ; 1D latitude OR 2D + yy, $ ; 1D longitude OR 2D +; 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_user, $ ; title of the plot ('' is default) + format=format ; format of colorbar annotations ('(F6.2)' is default) + + +PS_END, /PNG + +end Index: /trunk/mesoscale/PLOT/SPEC/SERIES/mawd_inc.pro =================================================================== --- /trunk/mesoscale/PLOT/SPEC/SERIES/mawd_inc.pro (revision 91) +++ /trunk/mesoscale/PLOT/SPEC/SERIES/mawd_inc.pro (revision 91) @@ -0,0 +1,35 @@ +;;;;;*************************************** PLOT TITLES + title_user = 'Water vapor column [precip. microns]' + title_axis = ['L!DS!N','North Latitude'] + +;;;;;*************************************** COLOR TABLES + flag_cb = 'true' + poscb = 0.75 + format = '(I0)' + colors = 32 + pal = 22;4;33 ;; GOOD: 4, 18, 22, 16, 37, 33, 39, 6, 11 + +;;;;;*************************************** FILL LIMITS + minfield_init = 0. + maxfield_init = 65. + ndiv = 13 + + ;;;;;*************************************** LIMIT TRICKS + ;;;;;********************* must always follow FILL LIMITS + missing_value=1.e30 + lim_max = maxfield_init & w=where((what_I_plot ge lim_max) and (what_I_plot le 1e9)) & if (w[0] ne -1) then what_I_plot[w]=lim_max + lim_min = minfield_init & w=where(what_I_plot le lim_min) & if (w[0] ne -1) then what_I_plot[w]=lim_min + ;lim_blank = 0.1 & w=where(abs(what_I_plot) le lim_blank) & if (w[0] ne -1) then what_I_plot[w]=missing_value + +;;;;;*************************************** CONTOUR + overcontour = overcontour + lev = minfield_init + (maxfield_init - minfield_init)*findgen(ndiv+1)/float(ndiv) + +;;;;;*************************************** AXIS + isotropic = 'true' + intervalx = 15.0 + intervaly = 15.0 + +;;;;;*************************************** MAP LIMITS + windowx = [30.,150.] + windowy = [-15.,90.] Index: unk/mesoscale/TESTS/newphys_tracers/callphys.def.mars11.radact =================================================================== --- /trunk/mesoscale/TESTS/newphys_tracers/callphys.def.mars11.radact (revision 90) +++ (revision ) @@ -1,116 +1,0 @@ -##General options -##~~~~~~~~~~~~~~~ -#Run with or without tracer transport ? -tracer=.true. - -#Diurnal cycle ? if diurnal=False, diurnal averaged solar heating -diurnal=.true. - -#Seasonal cycle ? if season=False, Ls stays constant, to value set in "start" -season = .true. - -#write some more output on the screen ? -lwrite = .false. - -#Save statistics in file "stats.nc" ? -callstats =.false. - -#Save EOF profiles in file "profiles" for Climate Database? -calleofdump = .false. - -## Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F) -## ~~~~~~~~~~~~~ -# =1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario, -# =4 Mars Year 24 from TES assimilation (old version of MY24; dust_tes.nc file) -# =24 Mars Year 24 from TES assimilation (ie: MCD reference case) -# =25 Mars Year 25 from TES assimilation (ie: a year with a global dust storm) -# =26 Mars Year 26 from TES assimilation -iaervar = 26 -# Dust opacity at 700 Pa (when constant, i.e. for the iaervar=1 case) -tauvis=0.2 -# Dust vertical distribution: -# (=0: old distrib. (Pollack90), =1: top set by "topdustref", -# =2: Viking scenario; =3 MGS scenario) -iddist = 3 -# Dust top altitude (km). (Matters only if iddist=1) -topdustref = 55. - -## Physical Parameterizations : -## ~~~~~~~~~~~~~~~~~~~~~~~~~~ -# call radiative transfer ? -callrad = .true. -# call NLTE radiative schemes ? matters only if callrad=T -callnlte = .false. -# call CO2 NIR absorption ? matters only if callrad=T -callnirco2 = .true. -# call turbulent vertical diffusion ? -calldifv = .true. -# call convective adjustment ? -calladj = .true. -# call CO2 condensation ? -callcond =.true. -# call thermal conduction in the soil ? -callsoil = .true. -# call Lott's gravity wave/subgrid topography scheme ? -calllott = .false. -# Impose polar cap surface albedos as observed by TES? -TESicealbedo = .true. - -## Radiative transfer options : -## ~~~~~~~~~~~~~~~~~~~~~~~~~~ -# the rad.transfer is computed every "iradia" physical timestep -iradia = 6 -# Output of the exchange coefficient mattrix ? for diagnostic only -callg2d = .false. -# Rayleigh scattering : (should be .false. for now) -rayleigh = .false. - -## Tracer (dust water, ice and/or chemical species) options (used if tracer=T): -## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# DUST: Transported dust ? (if >0, use 'dustbin' dust bins) -dustbin = 2 -# DUST: Radiatively active dust ? (matters if dustbin>0) -active = .true. -# DUST: use mass and number mixing ratios to predict dust size ? -doubleq = .true. -# DUST: lifted by GCM surface winds ? -lifting = .true. -# DUST: lifted by dust devils ? -callddevil = .false. -# DUST: Scavenging by CO2 snowfall ? -scavenging = .false. -# DUST/WATERICE: Gravitationnal sedimentation ? -sedimentation = .true. -# WATERICE: Radiatively active transported atmospheric water ice ? -activice = .true. -# WATER: Compute water cycle -water = .true. -# WATER: current permanent caps at both poles. True IS RECOMMENDED -# (with .true., North cap is a source of water and South pole -# is a cold trap) -caps = .false. #.true. -# PHOTOCHEMISTRY: include chemical species -photochem = .false. - -## Thermospheric options (relevant if tracer=T) : -##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# call thermosphere ? -callthermos = .false. -# WATER: included without cycle (only if water=.false.) -thermoswater = .false. -# call thermal conduction ? (only if callthermos=.true.) -callconduct = .false. -# call EUV heating ? (only if callthermos=.true.) -calleuv=.false. -# call molecular viscosity ? (only if callthermos=.true.) -callmolvis = .false. -# call molecular diffusion ? (only if callthermos=.true.) -callmoldiff = .false. -# call thermospheric photochemistry ? (only if callthermos=.true.) -thermochem = .false. -# date for solar flux calculation: (1985 < date < 2002) -## (Solar min=1996.4 ave=1993.4 max=1990.6) -solarcondate = 1993.4 - - -