source: trunk/MESOSCALE_DEV/PLOT/SPEC/LES/turbulence2.pro @ 211

Last change on this file since 211 was 85, checked in by aslmd, 14 years ago

LMD_MM_MARS et LMD_LES_MARS: ajout des routines IDL pour tracer les sorties --> voir mesoscale/PLOT
File size: 8.9 KB
Line 
1pro turbulence2
2
3
4retrieve='true'
5
6fast='true'
7;fast='false'
8
9
10p0=610. & t0=220. & r_cp=1/4.4 & grav=3.72 & R=192. & g_cp=1000.*grav/844.6
11;;;;;;;;;;;;;;;;;;;;;;;
12;                     ;
13; HINSON - HINSON !!! ;
14;                     ;
15r_cp = 0.25           ;
16;g_cp = 4.9           ;
17;;;;;;;;;;;;;;;;;;;;;;;
18
19
20
21;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
22;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
23;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
24;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
25
26
27
28;
29; input files
30;
31OPENR, 22, 'input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22
32OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23
33print, lonu, latu, lsu, lctu, hgtu, tsurfu
34
35;
36; retrieve
37;
38if (retrieve eq 'true') then begin
39
40;
41; get fields and dimensions
42;
43getcdf, file='t.nc', charvar='T', invar=t
44getcdf, file='p.nc', charvar='PTOT', invar=p
45getcdf, file='ph.nc', charvar='PHTOT', invar=ph
46nx=n_elements(t(*,0,0,0))
47ny=n_elements(t(0,*,0,0))
48nz=n_elements(t(0,0,*,0))
49nt=n_elements(t(0,0,0,*))
50   ;;;; pour eviter de detecter des fluctuations plus basses que PBL top ;;;
51     smooth_gridpoint = 2 & smooth_time = 2 & smooth_vert = 2 
52     ;smooth_gridpoint = 0 & smooth_time = 0 & smooth_vert = 0 
53     t = smooth(t,[smooth_gridpoint,smooth_gridpoint,smooth_vert,smooth_time],/EDGE_TRUNCATE)
54     p = smooth(p,[smooth_gridpoint,smooth_gridpoint,smooth_vert,smooth_time],/EDGE_TRUNCATE)
55     ;ph = smooth(ph,[smooth_gridpoint,smooth_gridpoint,smooth_vert,smooth_time],/EDGE_TRUNCATE)
56   ;;;; pour eviter de detecter des fluctuations plus basses que PBL top ;;;
57
58;
59; localtime loop
60;
61localtime = lctu + 100.*findgen(nt)/3700.
62  ;
63  ; radio-occultations
64  ;
65  user_lt=17.0 & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
66  user_lt=17.5 & yeah2=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
67  ;;
68  ;; PBL growth
69  ;;
70  ;user_lt=14.5 & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
71  ;user_lt=18.0 & yeah2=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
72  ;
73
74  user_lt=14.0 & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
75  user_lt=14.5 & yeah2=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt))))
76
77
78
79deb=yeah(0)
80fin=yeah2(0)
81
82;
83; calculate over domain
84;
85pbl_depth_tab=findgen(nx,ny,fin-deb+1)
86nloop = nx*ny
87
88;
89; LOOOOOOP
90;
91for tt=deb,fin do begin
92
93 tprof_mean = 0.
94 heightp_mean = 0.
95 staticstab_mean = 0.
96 nbad = 0
97
98for i=0,nx-1 do begin
99for j=0,ny-1 do begin
100
101 height=reform( ph( i,j,*,tt ) ) / 1000. / grav
102 heightp=height(0:nz-1) ;no stagger
103
104 thprof = t0 + reform( t( i,j,*,tt ) )
105 pprof = reform( p( i,j,*,tt ) )
106 tprof = thprof * ( pprof/p0 )^r_cp
107
108 staticstab = DERIV( heightp , tprof ) + g_cp 
109
110if (fast eq 'true') then begin
111
112; MOINS BON mais PLUS RAPIDE
113tprof_mean = tprof_mean + tprof/nloop
114heightp_mean = heightp_mean + heightp/nloop
115staticstab_mean = staticstab_mean + staticstab/nloop
116
117endif else begin
118
119;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
120;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
121;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
122  z_in = heightp
123  profile_in = staticstab
124  resol = 500.  ;200 un peu juste
125  ;********************************************************************
126  ; ---- numerical recipes in C - section 3.3 ----
127  ;
128  ; Calculate interpolating cubic spline
129        yspline = SPL_INIT(z_in,profile_in)
130  ; Define the X values P at which we desire interpolated Y values
131        xspline=min(z_in)+findgen(floor(max(z_in)-min(z_in))*resol)/resol
132  ; Calculate the interpolated Y values
133        result=spl_interp(z_in,profile_in,yspline,xspline)
134  ;********************************************************************
135  staticstab = result
136  heightp = xspline
137
138  heightp_mean = heightp_mean + heightp/nloop
139  staticstab_mean = staticstab_mean + staticstab/nloop
140    yspline = SPL_INIT(z_in,tprof)
141    tprof = spl_interp(z_in,tprof,yspline,xspline)
142  tprof_mean = tprof_mean + tprof/nloop
143;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
144;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
145;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
146
147endelse
148
149lim=3.
150lim=1.5
151w = where( (staticstab ge lim) and (heightp-heightp(0) ge lim) )
152pbl_depth_tab(i,j,tt-deb) = heightp(w(0)) - heightp(0)
153
154   ; garde-fou supplementaire
155   if (heightp(w(0)) - heightp(0) lt 2.) then begin
156     pbl_depth_tab(i,j,tt-deb) = !VALUES.F_NAN
157     ;plot, profile_in, z_in, xrange=[-1.,5.], psym=5, title='bad!' & oplot, staticstab, heightp & oplot, staticstab, findgen(nz)*0. + heightp(w(0)), linestyle=2
158     print, heightp(w) - heightp(0)
159     nbad = nbad+1
160   endif
161
162;plot, profile_in, z_in, xrange=[-1.,5.], psym=5 & oplot, staticstab, heightp & oplot, staticstab, findgen(nz)*0. + heightp(w(0)), linestyle=2
163;pause
164
165endfor
166;print, i
167endfor
168
169;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
170;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
171  print, '-------------------------------------------------------'
172  print, 'Local Time : ', localtime(tt)
173  print, 'nbad ', nbad
174  print, 'PBL height - min ', min(pbl_depth_tab(*,*,tt-deb))
175  print, 'PBL height - max ', max(pbl_depth_tab(*,*,tt-deb))
176  print, 'PBL height - mean ', mean(pbl_depth_tab(*,*,tt-deb))
177    ;; static stability of mean profile
178    staticstab = DERIV( heightp_mean , tprof_mean ) + g_cp
179
180          ;z_in = heightp_mean
181          ;profile_in = staticstab
182          ;resol = 500  ;200 un peu juste
183          ;;********************************************************************
184          ;; ---- numerical recipes in C - section 3.3 ----
185          ;;
186          ;; Calculate interpolating cubic spline
187          ;      yspline = SPL_INIT(z_in,profile_in)
188          ;; Define the X values P at which we desire interpolated Y values
189          ;      xspline=min(z_in)+findgen(floor(max(z_in))*resol)/resol
190          ;; Calculate the interpolated Y values
191          ;      result=spl_interp(z_in,profile_in,yspline,xspline)
192          ;;********************************************************************
193          ;staticstab = result
194          ;w = where((staticstab ge 1.5) and (xspline-xspline(0) gt 2.))
195          ;print, 'PBL height - mean profile', xspline(w(0)) - xspline(0)
196
197    w = where((staticstab ge 1.5) and (heightp_mean-heightp_mean(0) gt 1.5))
198    print, 'PBL height - mean profile', heightp_mean(w(0))-heightp_mean(0)
199    ;; mean static stability
200    w = where((staticstab_mean ge 1.5) and (heightp_mean-heightp_mean(0) gt 1.5))
201    print, 'PBL height - mean static stability', heightp_mean(w(0)) - heightp_mean(0)
202
203   ;plot, staticstab, heightp_mean, xrange=[-1.,5.]
204   ;oplot, staticstab, findgen(nz)*0. + heightp_mean(w(0)), linestyle=2
205   ;oplot, staticstab_mean, heightp_mean
206   ;oplot, staticstab_mean, findgen(nz)*0. + heightp_mean(w(0)), linestyle=2
207   ;pause
208;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
209;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
210endfor
211endif
212
213  print, '-------------------------------------------------------'
214  print, 'PBL height - min ', min(pbl_depth_tab(*,*,*))
215  print, 'PBL height - max ', max(pbl_depth_tab(*,*,*))
216  print, 'PBL height - mean ', mean(pbl_depth_tab(*,*,*))
217  print, '-------------------------------------------------------'
218
219stop
220
221!p.charthick = 2.0
222!p.thick = 3.0
223!x.thick = 2.0
224!y.thick = 2.0
225set_plot, 'ps' & device, filename='plot/map.ps', /color
226;map_latlon, reform(pbl_depth_tab(*,*,0))
227
228
229what_I_plot = reform(pbl_depth_tab(*,*,0))
230pal=4
231colors=128
232title_user='PBL depth'
233map_latlon, $
234        what_I_plot, $                          ; 2D field
235        lon, $                                  ; 1D latitude
236        lat, $                                  ; 1D longitude
237;        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
238;        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
239;        overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
240;        overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
241;        overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
242        ct=pal, $                               ; color table (33-rainbow is default)
243        colors=colors, $                        ; number of colors/levels (32 is default)
244        title=title_user;, $                     ; title of the plot ('' is default)
245;        format=format                           ; format of colorbar annotations ('(F6.2)' is default)
246
247
248
249device, /close
250end
Note: See TracBrowser for help on using the repository browser.