| [85] | 1 | pro turbulence2 |
|---|
| 2 | |
|---|
| 3 | |
|---|
| 4 | retrieve='true' |
|---|
| 5 | |
|---|
| 6 | fast='true' |
|---|
| 7 | ;fast='false' |
|---|
| 8 | |
|---|
| 9 | |
|---|
| 10 | p0=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 | ; ; |
|---|
| 15 | r_cp = 0.25 ; |
|---|
| 16 | ;g_cp = 4.9 ; |
|---|
| 17 | ;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 18 | |
|---|
| 19 | |
|---|
| 20 | |
|---|
| 21 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 22 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 23 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 24 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 25 | |
|---|
| 26 | |
|---|
| 27 | |
|---|
| 28 | ; |
|---|
| 29 | ; input files |
|---|
| 30 | ; |
|---|
| 31 | OPENR, 22, 'input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22 |
|---|
| 32 | OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23 |
|---|
| 33 | print, lonu, latu, lsu, lctu, hgtu, tsurfu |
|---|
| 34 | |
|---|
| 35 | ; |
|---|
| 36 | ; retrieve |
|---|
| 37 | ; |
|---|
| 38 | if (retrieve eq 'true') then begin |
|---|
| 39 | |
|---|
| 40 | ; |
|---|
| 41 | ; get fields and dimensions |
|---|
| 42 | ; |
|---|
| 43 | getcdf, file='t.nc', charvar='T', invar=t |
|---|
| 44 | getcdf, file='p.nc', charvar='PTOT', invar=p |
|---|
| 45 | getcdf, file='ph.nc', charvar='PHTOT', invar=ph |
|---|
| 46 | nx=n_elements(t(*,0,0,0)) |
|---|
| 47 | ny=n_elements(t(0,*,0,0)) |
|---|
| 48 | nz=n_elements(t(0,0,*,0)) |
|---|
| 49 | nt=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 | ; |
|---|
| 61 | localtime = 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 | |
|---|
| 79 | deb=yeah(0) |
|---|
| 80 | fin=yeah2(0) |
|---|
| 81 | |
|---|
| 82 | ; |
|---|
| 83 | ; calculate over domain |
|---|
| 84 | ; |
|---|
| 85 | pbl_depth_tab=findgen(nx,ny,fin-deb+1) |
|---|
| 86 | nloop = nx*ny |
|---|
| 87 | |
|---|
| 88 | ; |
|---|
| 89 | ; LOOOOOOP |
|---|
| 90 | ; |
|---|
| 91 | for tt=deb,fin do begin |
|---|
| 92 | |
|---|
| 93 | tprof_mean = 0. |
|---|
| 94 | heightp_mean = 0. |
|---|
| 95 | staticstab_mean = 0. |
|---|
| 96 | nbad = 0 |
|---|
| 97 | |
|---|
| 98 | for i=0,nx-1 do begin |
|---|
| 99 | for 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 | |
|---|
| 110 | if (fast eq 'true') then begin |
|---|
| 111 | |
|---|
| 112 | ; MOINS BON mais PLUS RAPIDE |
|---|
| 113 | tprof_mean = tprof_mean + tprof/nloop |
|---|
| 114 | heightp_mean = heightp_mean + heightp/nloop |
|---|
| 115 | staticstab_mean = staticstab_mean + staticstab/nloop |
|---|
| 116 | |
|---|
| 117 | endif 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 | |
|---|
| 147 | endelse |
|---|
| 148 | |
|---|
| 149 | lim=3. |
|---|
| 150 | lim=1.5 |
|---|
| 151 | w = where( (staticstab ge lim) and (heightp-heightp(0) ge lim) ) |
|---|
| 152 | pbl_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 | |
|---|
| 165 | endfor |
|---|
| 166 | ;print, i |
|---|
| 167 | endfor |
|---|
| 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 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 210 | endfor |
|---|
| 211 | endif |
|---|
| 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 | |
|---|
| 219 | stop |
|---|
| 220 | |
|---|
| 221 | !p.charthick = 2.0 |
|---|
| 222 | !p.thick = 3.0 |
|---|
| 223 | !x.thick = 2.0 |
|---|
| 224 | !y.thick = 2.0 |
|---|
| 225 | set_plot, 'ps' & device, filename='plot/map.ps', /color |
|---|
| 226 | ;map_latlon, reform(pbl_depth_tab(*,*,0)) |
|---|
| 227 | |
|---|
| 228 | |
|---|
| 229 | what_I_plot = reform(pbl_depth_tab(*,*,0)) |
|---|
| 230 | pal=4 |
|---|
| 231 | colors=128 |
|---|
| 232 | title_user='PBL depth' |
|---|
| 233 | map_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 | |
|---|
| 249 | device, /close |
|---|
| 250 | end |
|---|
