| [85] | 1 | pro getturb, saveps=saveps, noplot=noplot |
|---|
| 2 | |
|---|
| 3 | ;---------------------------------- |
|---|
| 4 | ; USE: getturb |
|---|
| 5 | ; getturb, saveps='false' |
|---|
| 6 | ;---------------------------------- |
|---|
| 7 | |
|---|
| 8 | history_interval_s = 400. |
|---|
| 9 | history_interval_s = 100. |
|---|
| 10 | smoothampl=3700/history_interval_s |
|---|
| 11 | ;smoothampl=0. ;; no smooth |
|---|
| 12 | ;smoothampl=10. |
|---|
| 13 | smoothampl=2. |
|---|
| 14 | ;smoothampl=5. |
|---|
| 15 | |
|---|
| 16 | ; |
|---|
| 17 | ; constantes |
|---|
| 18 | ; |
|---|
| 19 | p0=610. & t0=220. & r_cp=1/4.4 & grav=3.72 & R=192. |
|---|
| 20 | |
|---|
| 21 | ; |
|---|
| 22 | ; graphics definition |
|---|
| 23 | ; |
|---|
| 24 | if (n_elements(saveps) eq 0) then saveps='true' |
|---|
| 25 | if (n_elements(noplot) eq 0) then noplot='false' |
|---|
| 26 | if (saveps eq 'false') then begin |
|---|
| 27 | ;!p.multi=[0,3,2] |
|---|
| 28 | !P.CHARSIZE=2. |
|---|
| 29 | WINDOW, /PIXMAP & WDELETE & DEVICE,BYPASS_TRANSLATION=0,DECOMPOSED=0,RETAIN=2 |
|---|
| 30 | endif else begin |
|---|
| 31 | PREF_SET, 'IDL_PATH', '/home/spiga/Save/SOURCES/IDL/fsc_psconfig:<IDL_DEFAULT>', /COMMIT |
|---|
| 32 | endelse |
|---|
| 33 | |
|---|
| 34 | ; |
|---|
| 35 | ; retrieve fields |
|---|
| 36 | ; |
|---|
| 37 | openr,unit,'getturb.dat',/get_lun,error=err |
|---|
| 38 | IF (err ne 0) THEN BEGIN |
|---|
| 39 | |
|---|
| 40 | ; |
|---|
| 41 | ; input files |
|---|
| 42 | ; |
|---|
| 43 | OPENR, 22, 'input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22 |
|---|
| 44 | OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23 |
|---|
| 45 | |
|---|
| 46 | ; |
|---|
| 47 | ; get fields |
|---|
| 48 | ; |
|---|
| 49 | domain='d01' & filesWRF = FindFile('wrfout_'+domain+'_????-??-??_??:??:??') & nf=n_elements(filesWRF) |
|---|
| 50 | |
|---|
| 51 | |
|---|
| 52 | ; get dimensions |
|---|
| 53 | ; |
|---|
| 54 | id=ncdf_open(filesWRF(0)) |
|---|
| 55 | NCDF_DIMINQ, id, NCDF_DIMID(id, 'west_east' ), toto, nx & NCDF_DIMINQ, id, NCDF_DIMID(id, 'south_north' ), toto, ny |
|---|
| 56 | NCDF_DIMINQ, id, NCDF_DIMID(id, 'bottom_top' ), toto, nz & NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), toto, nt |
|---|
| 57 | NCDF_CLOSE, id |
|---|
| 58 | id=ncdf_open(filesWRF(nf-1)) ;; for interrupted runs |
|---|
| 59 | NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), toto, ntlast |
|---|
| 60 | NCDF_CLOSE, id |
|---|
| 61 | |
|---|
| 62 | ; |
|---|
| 63 | ; prepare loop |
|---|
| 64 | ; |
|---|
| 65 | ;nloop1 = nf-1 & nloop2 = nt & yeye = 0 |
|---|
| 66 | yeye = 0 & nttot = (nf-1)*nt + ntlast |
|---|
| 67 | |
|---|
| 68 | localtime = lctu + history_interval_s*findgen(nttot)/3700. |
|---|
| 69 | |
|---|
| 70 | wt = fltarr(nz,nttot) & tke = fltarr(nz,nttot) & ztke = fltarr(nz,nttot) & t = fltarr(nz,nttot) |
|---|
| 71 | p = fltarr(nz) & ph = fltarr(nz) & pht = fltarr(nz,nttot) & pt = fltarr(nz,nttot) & stst = fltarr(nz,nttot) |
|---|
| 72 | xtke = fltarr(nz,nttot) & ytke = fltarr(nz,nttot) & wmax = fltarr(nz,nttot) & wmin = fltarr(nz,nttot) |
|---|
| 73 | depressions = fltarr(nttot) & psmin = fltarr(nttot) |
|---|
| 74 | |
|---|
| 75 | ; |
|---|
| 76 | ; loop loop |
|---|
| 77 | ; |
|---|
| 78 | for loop = 0, nf-1 do begin |
|---|
| 79 | timetime = SYSTIME(1) |
|---|
| 80 | if (loop ne nf-1) then nloop2=nt else nloop2=ntlast |
|---|
| 81 | |
|---|
| 82 | for loop2 = 0, nloop2-1 do begin |
|---|
| 83 | |
|---|
| 84 | anomalt = 1. & anomalu = 1. & anomalv = 1. |
|---|
| 85 | ;print, loop, loop2 |
|---|
| 86 | |
|---|
| 87 | ; ------------ |
|---|
| 88 | ; t' = t - <t> |
|---|
| 89 | ; ------------ |
|---|
| 90 | tprime = getget(filesWRF(loop), 'T', anomaly=anomalt, count=[0,0,0,1], offset=[0,0,0,loop2]) ;; t' = t - <t> |
|---|
| 91 | t(*,yeye) = t0 + TEMPORARY(anomalt) |
|---|
| 92 | ; ------ |
|---|
| 93 | ; w' = w |
|---|
| 94 | ; ------ |
|---|
| 95 | wprime = getget(filesWRF(loop), 'W', count=[0,0,0,1], offset=[0,0,0,loop2]) |
|---|
| 96 | for toto = 0, nz-1 do begin |
|---|
| 97 | wmax(toto,yeye) = max(reform(wprime(*,*,toto,0)),min=tutu) & wmin(toto,yeye) = tutu |
|---|
| 98 | endfor |
|---|
| 99 | ; ; ------ |
|---|
| 100 | ; ; u' = u and v' = v (car PAS de background wind !) |
|---|
| 101 | ; ; ------ |
|---|
| 102 | ; uprime = getget(filesWRF(loop), 'U', count=[0,0,0,1], offset=[0,0,0,loop2]) |
|---|
| 103 | ; vprime = getget(filesWRF(loop), 'V', count=[0,0,0,1], offset=[0,0,0,loop2]) |
|---|
| 104 | ; -------------------------------------------------------- |
|---|
| 105 | ; tke = 0.5 ( <u'^2> + <v'^2> + <w'^2> ) ; u' = u ; v' = v |
|---|
| 106 | ; -------------------------------------------------------- |
|---|
| 107 | ; ztke is 0.5 * <w'^2>/2 or 0.5 * sigma_w^2 |
|---|
| 108 | ztke(*,yeye) = 0.5 * TOTAL(TOTAL(wprime^2,1),1) / float(nx) / float(ny) |
|---|
| 109 | xtke(*,yeye) = 0.5 * TOTAL(TOTAL(getget(filesWRF(loop), 'U', anomaly=anomalu, count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) / float(nx) / float(ny) |
|---|
| 110 | ytke(*,yeye) = 0.5 * TOTAL(TOTAL(getget(filesWRF(loop), 'V', anomaly=anomalv, count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) / float(nx) / float(ny) |
|---|
| 111 | ; xtke(*,yeye) = 0.5 * TOTAL(TOTAL(uprime^2,1),1) / float(nx) / float(ny) |
|---|
| 112 | ; ytke(*,yeye) = 0.5 * TOTAL(TOTAL(vprime^2,1),1) / float(nx) / float(ny) |
|---|
| 113 | tke(*,yeye) = xtke(*,yeye) + $ |
|---|
| 114 | ytke(*,yeye) + $ |
|---|
| 115 | ztke(*,yeye) |
|---|
| 116 | ; ------ |
|---|
| 117 | ; <w't'> |
|---|
| 118 | ; ------ |
|---|
| 119 | wt(*,yeye) = TOTAL(TOTAL(TEMPORARY(tprime) * TEMPORARY(wprime),1),1) / float(nx) / float(ny) |
|---|
| 120 | ; ------ |
|---|
| 121 | ; p & ph |
|---|
| 122 | ; ------ |
|---|
| 123 | ;if (loop + loop2 eq 0) then nloopbeware = nloop2-1 else nloopbeware = nloop2 ; 1ere valeur vaut 0 |
|---|
| 124 | nttotbeware = nttot-1 ; 1ere valeur vaut 0 |
|---|
| 125 | pht(*,yeye) = TOTAL(TOTAL(getget(filesWRF(loop), 'PHTOT', count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny) / 1000. / 3.72 |
|---|
| 126 | pt(*,yeye) = TOTAL(TOTAL(getget(filesWRF(loop), 'PTOT' , count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny) |
|---|
| 127 | ph = TEMPORARY(ph) + pht(*,yeye) / nttotbeware ;/ nloopbeware / nloop1 |
|---|
| 128 | p = TEMPORARY(p ) + pt(*,yeye) / nttot ;/ nloop2 / nloop1 |
|---|
| 129 | ; ---------------- |
|---|
| 130 | ; static stability |
|---|
| 131 | ; ---------------- |
|---|
| 132 | stst(*,yeye) = DERIV( reform(pht(*,yeye)) - hgtu/1000. , reform(t(*,yeye)) * ( reform(pt(*,yeye)) /p0 )^r_cp ) + 1000.*grav / (R / r_cp) ;; 4.9 dans Hinson |
|---|
| 133 | ; ---------------- |
|---|
| 134 | ; surface pressure |
|---|
| 135 | ; ---------------- |
|---|
| 136 | !QUIET=1 |
|---|
| 137 | ;psfc = getget(filesWRF(loop), 'PSFC', anomaly=1, count=[0,0,0,1], offset=[0,0,0,loop2]) ;; plus rapide que autre routine! pas grave avertissement |
|---|
| 138 | psfc = getget(filesWRF(loop), 'PSFC', count=[0,0,0,1], offset=[0,0,0,loop2]) |
|---|
| 139 | !QUIET=0 |
|---|
| 140 | ; psfc = getget2d(filesWRF(loop), 'PSFC', anomaly=1, count=[0,0,1], offset=[0,0,loop2]) |
|---|
| 141 | psmin(yeye) = min(psfc) ;& print, psmin(yeye) |
|---|
| 142 | w = where(TEMPORARY(psfc) lt -0.5) ;& print, n_elements(w) |
|---|
| 143 | if (w(0) ne -1) then depressions(yeye) = 1000. * n_elements(w) / float(nx) / float(ny) else depressions(yeye) = 0. |
|---|
| 144 | ; ---------- |
|---|
| 145 | ; loop count |
|---|
| 146 | ; ---------- |
|---|
| 147 | yeye = TEMPORARY(yeye) + 1 ;& print, yeye & print, SYSTIME(1) - timetime, ' s' |
|---|
| 148 | |
|---|
| 149 | endfor |
|---|
| 150 | print, 'file '+string(loop+1,'(I0)'), SYSTIME(1) - timetime, ' s' |
|---|
| 151 | endfor |
|---|
| 152 | h = TEMPORARY(ph) - hgtu/1000. ;; altitude above ground |
|---|
| 153 | ht = TEMPORARY(pht) - hgtu/1000. |
|---|
| 154 | ; |
|---|
| 155 | ; save |
|---|
| 156 | ; |
|---|
| 157 | save, wt, tke, ztke, h, ht, t, p, pt, stst, localtime, xtke, ytke, wmax, wmin, depressions, psmin, filename='getturb.dat' |
|---|
| 158 | nz = n_elements(h) |
|---|
| 159 | nnn = n_elements(tke(0,*)) |
|---|
| 160 | |
|---|
| 161 | ENDIF ELSE BEGIN |
|---|
| 162 | |
|---|
| 163 | print, 'OK, file is here' |
|---|
| 164 | restore, filename='getturb.dat' |
|---|
| 165 | nz = n_elements(h) |
|---|
| 166 | nnn = n_elements(tke(0,*)) |
|---|
| 167 | |
|---|
| 168 | OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23 |
|---|
| 169 | |
|---|
| 170 | @include_ustar.pro |
|---|
| 171 | |
|---|
| 172 | ;@include_hfx.pro |
|---|
| 173 | |
|---|
| 174 | ;@include_w.pro |
|---|
| 175 | |
|---|
| 176 | ENDELSE |
|---|
| 177 | |
|---|
| 178 | ; |
|---|
| 179 | ; pbl height |
|---|
| 180 | ; |
|---|
| 181 | pbl_height = fltarr(nnn) & hfsurf = fltarr(nnn) & w_star = fltarr(nnn) & hfbot = fltarr(nnn) & w_star_bot = fltarr(nnn) & hfbotmax = fltarr(nnn) |
|---|
| 182 | a_vel = fltarr(nz,nnn) & a_h = fltarr(nz,nnn) & a_wt = fltarr(nz,nnn) & a_vel_bot = fltarr(nz,nnn) & a_wt_bot = fltarr(nz,nnn) |
|---|
| 183 | tt_top = fltarr(nnn) & z_top = fltarr(nnn) & t_bot = fltarr(nnn) & p_bot = fltarr(nnn) |
|---|
| 184 | tt_bot = fltarr(nnn) |
|---|
| 185 | wmaxmax = fltarr(nnn) |
|---|
| 186 | wminmin = fltarr(nnn) |
|---|
| 187 | tkemax = fltarr(nnn) |
|---|
| 188 | |
|---|
| 189 | for i=1, nnn-1 do begin ;; ne pas traiter l'init |
|---|
| 190 | |
|---|
| 191 | ; |
|---|
| 192 | ; pbl height |
|---|
| 193 | ; |
|---|
| 194 | z_in = reform(ht(*,i)) |
|---|
| 195 | profile_in = reform(stst(*,i)) ;; wt pas suffisamment regulier |
|---|
| 196 | resol = 1000 ;200 un peu juste |
|---|
| 197 | ;******************************************************************** |
|---|
| 198 | ; ---- numerical recipes in C - section 3.3 ---- |
|---|
| 199 | ; |
|---|
| 200 | ; Calculate interpolating cubic spline |
|---|
| 201 | yspline = SPL_INIT(z_in,profile_in) |
|---|
| 202 | ; Define the X values P at which we desire interpolated Y values |
|---|
| 203 | xspline=min(z_in)+findgen(floor(max(z_in))*resol)/resol |
|---|
| 204 | ; Calculate the interpolated Y values |
|---|
| 205 | result=spl_interp(z_in,profile_in,yspline,xspline) |
|---|
| 206 | ;******************************************************************** |
|---|
| 207 | w = where( (result ge 1.5) and (xspline gt 0.01) ) ;0.3) ) ;0.1 suffit si HR |
|---|
| 208 | ;if ( localtime(i) gt 16. ) then w = where( (result ge 1.5) and (xspline gt 2.) ) |
|---|
| 209 | |
|---|
| 210 | if ( localtime(i) gt 15. ) then w = where( (result ge 1.5) and (xspline gt 1.) ) |
|---|
| 211 | |
|---|
| 212 | ;;;; special tau = 5. |
|---|
| 213 | ;if ( localtime(i) ge 13. ) then w = where( (result ge 1.5) and (xspline gt 0.35) ) |
|---|
| 214 | ;if ( localtime(i) ge 16.05 ) then w = [0] |
|---|
| 215 | |
|---|
| 216 | pbl_height(i) = xspline(w(0)) ;& print, 'PBL height - mean profile', pbl_height(i) |
|---|
| 217 | z_top(i) = xspline(w(0)) + hgtu/1000. |
|---|
| 218 | |
|---|
| 219 | ; |
|---|
| 220 | ; temperature @ pbl height |
|---|
| 221 | ; |
|---|
| 222 | z_in = reform(ht(*,i)) |
|---|
| 223 | profile_in = reform(t(*,i)) |
|---|
| 224 | ;profile_in = reform(t(*,i))*(reform(pt(*,i))/p0)^r_cp |
|---|
| 225 | profile_in2 = reform(pt(*,i)) |
|---|
| 226 | resol = 1000 ;200 un peu juste |
|---|
| 227 | ;******************************************************************** |
|---|
| 228 | ; ---- numerical recipes in C - section 3.3 ---- |
|---|
| 229 | ; |
|---|
| 230 | ; Calculate interpolating cubic spline |
|---|
| 231 | yspline = SPL_INIT(z_in,profile_in) |
|---|
| 232 | yspline2 = SPL_INIT(z_in,profile_in2) |
|---|
| 233 | ; Define the X values P at which we desire interpolated Y values |
|---|
| 234 | xspline=min(z_in)+findgen(floor(max(z_in))*resol)/resol |
|---|
| 235 | ; Calculate the interpolated Y values |
|---|
| 236 | result=spl_interp(z_in,profile_in,yspline,xspline) |
|---|
| 237 | result2=spl_interp(z_in,profile_in2,yspline2,xspline) |
|---|
| 238 | ;******************************************************************** |
|---|
| 239 | w = where( abs(xspline-pbl_height(i)) eq min(abs(xspline-pbl_height(i))) ) & tt_top(i) = result[w(0)]*(result2[w(0)]/p0)^r_cp ;& print, z_top(i), t_top(i) |
|---|
| 240 | hhh = 1000. / 1000. ;100. / 1000. |
|---|
| 241 | w = where( abs(xspline-hhh) eq min(abs(xspline-hhh)) ) & t_bot(i) = result[w(0)] & p_bot(i) = result2[w(0)] ;& print, p_bot(i), t_bot(i) |
|---|
| 242 | tt_bot(i) = result[w(0)]*(result2[w(0)]/p0)^r_cp |
|---|
| 243 | |
|---|
| 244 | ; |
|---|
| 245 | ; heat flux at the "surface" |
|---|
| 246 | ; |
|---|
| 247 | hfsurf(i) = wt(1,i) ;; attention vaut 0 a la hauteur 0 |
|---|
| 248 | |
|---|
| 249 | ; |
|---|
| 250 | ; heat flux at the bottom of mixed layer = top of radiative layer |
|---|
| 251 | ; |
|---|
| 252 | z_in = reform(ht(*,i)) |
|---|
| 253 | profile_in = reform(wt(*,i)) ;; wt pas suffisamment regulier |
|---|
| 254 | resol = 1000 ;200 un peu juste |
|---|
| 255 | ;******************************************************************** |
|---|
| 256 | ; ---- numerical recipes in C - section 3.3 ---- |
|---|
| 257 | ; |
|---|
| 258 | ; Calculate interpolating cubic spline |
|---|
| 259 | yspline = SPL_INIT(z_in,profile_in) |
|---|
| 260 | ; Define the X values P at which we desire interpolated Y values |
|---|
| 261 | xspline=min(z_in)+findgen(floor(max(z_in))*resol)/resol |
|---|
| 262 | ; Calculate the interpolated Y values |
|---|
| 263 | result=spl_interp(z_in,profile_in,yspline,xspline) |
|---|
| 264 | ;******************************************************************** |
|---|
| 265 | hhh = 300./1000. & yeah=where(abs(xspline-hhh) eq (min(abs(xspline-hhh)))) & nnzz=yeah(0) & hfbot(i)=result(nnzz) ;; 300m semble OK, 500m super |
|---|
| 266 | w=where(xspline lt 1.5) & result=result[w] & xspline=xspline[w] & yeah=where(result eq max(result)) & nnzz=yeah(0) & hfbotmax(i)=result(nnzz) |
|---|
| 267 | ;print, 'surf, 300m, max ', hfsurf(i), hfbot(i), hfbotmax(i) ;& plot, result, xspline |
|---|
| 268 | |
|---|
| 269 | result = reform(wmax(*,i)) & yeah=where(result eq max(result)) & nnzz=yeah(0) & wmaxmax(i)=result(nnzz) |
|---|
| 270 | result = reform(wmin(*,i)) & yeah=where(result eq min(result)) & nnzz=yeah(0) & wminmin(i)=result(nnzz) |
|---|
| 271 | result = reform(tke(*,i)) & yeah=where(result eq max(result)) & nnzz=yeah(0) & tkemax(i )=result(nnzz) |
|---|
| 272 | |
|---|
| 273 | ; |
|---|
| 274 | ; free convection velocity scale |
|---|
| 275 | ; |
|---|
| 276 | ; w_star_bot(i) = ( 1000.* grav * pbl_height(i) * hfbot(i) / t(0,i) ) ^ (1./3.) |
|---|
| 277 | w_star_bot(i) = ( 1000.* grav * pbl_height(i) * hfbotmax(i) / t(0,i) ) ^ (1./3.) |
|---|
| 278 | w_star(i) = ( 1000.* grav * pbl_height(i) * hfsurf(i) / t(0,i) ) ^ (1./3.) |
|---|
| 279 | |
|---|
| 280 | ; |
|---|
| 281 | ; dimensional quantities |
|---|
| 282 | ; |
|---|
| 283 | a_vel(*,i) = reform( 2.*ztke(*,i) / w_star(i)^2 ) |
|---|
| 284 | a_vel_bot(*,i) = reform( 2.*ztke(*,i) / w_star_bot(i)^2 ) |
|---|
| 285 | a_h(*,i) = reform( ht(*,i) / pbl_height(i) ) |
|---|
| 286 | a_wt(*,i) = reform( wt(*,i) / hfsurf(i) ) ;reform( wt(*,i) / w_star(i) ) |
|---|
| 287 | a_wt_bot(*,i) = reform( wt(*,i) / hfbotmax(i) ) ;reform( wt(*,i) / w_star_bot(i) ) |
|---|
| 288 | |
|---|
| 289 | endfor |
|---|
| 290 | ; |
|---|
| 291 | ; free convection time scale |
|---|
| 292 | ; |
|---|
| 293 | fcts = 1000. * pbl_height / w_star / 3700. |
|---|
| 294 | ; |
|---|
| 295 | ; mixed layer temperature scale |
|---|
| 296 | ; |
|---|
| 297 | mlts = hfsurf / w_star |
|---|
| 298 | |
|---|
| 299 | ;pbl_height = SMOOTH(TEMPORARY(pbl_height), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 300 | ;hfsurf = SMOOTH(TEMPORARY(hfsurf), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 301 | ;w_star = SMOOTH(TEMPORARY(w_star), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 302 | ;hfbot = SMOOTH(TEMPORARY(hfbot), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 303 | ;w_star_bot = SMOOTH(TEMPORARY(w_star_bot), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 304 | ;hfbotmax = SMOOTH(TEMPORARY(hfbotmax), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 305 | ;a_vel = SMOOTH(TEMPORARY(a_vel), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 306 | ;a_h = SMOOTH(TEMPORARY(a_h), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 307 | ;a_wt = SMOOTH(TEMPORARY(a_wt), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 308 | ;a_vel_bot = SMOOTH(TEMPORARY(a_vel_bot), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 309 | ;a_wt_bot = SMOOTH(TEMPORARY(a_wt_bot), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 310 | ;tt_top = SMOOTH(TEMPORARY(tt_top), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 311 | ;z_top = SMOOTH(TEMPORARY(z_top), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 312 | ;t_bot = SMOOTH(TEMPORARY(t_bot), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 313 | ;p_bot = SMOOTH(TEMPORARY(p_bot), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 314 | |
|---|
| 315 | ; |
|---|
| 316 | ; save for comparison sake |
|---|
| 317 | ; |
|---|
| 318 | save, pbl_height, hfbot, hfsurf, w_star, w_star_bot, a_vel_bot, a_vel, a_h, a_wt_bot, a_wt, fcts, mlts, hfbotmax, tt_top, z_top, t_bot, p_bot, tt_bot, wmaxmax, wminmin, tkemax, filename='compturb.dat' |
|---|
| 319 | if (noplot ne 'false') then stop |
|---|
| 320 | |
|---|
| 321 | ; |
|---|
| 322 | ; smooth smooth |
|---|
| 323 | ; |
|---|
| 324 | wt = SMOOTH(TEMPORARY(wt), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 325 | tke = SMOOTH(TEMPORARY(tke), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 326 | ztke = SMOOTH(TEMPORARY(ztke), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 327 | xtke = SMOOTH(TEMPORARY(xtke), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 328 | ytke = SMOOTH(TEMPORARY(ytke), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 329 | wmax = SMOOTH(TEMPORARY(wmax), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 330 | wmin = SMOOTH(TEMPORARY(wmin), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 331 | depressions = SMOOTH(TEMPORARY(depressions), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 332 | psmin = SMOOTH(TEMPORARY(psmin), [smoothampl], /EDGE_TRUNCATE) |
|---|
| 333 | |
|---|
| 334 | |
|---|
| 335 | ;*******************; |
|---|
| 336 | ;*******************; |
|---|
| 337 | ; PLOTS PLOTS PLOTS ; |
|---|
| 338 | ;*******************; |
|---|
| 339 | ;*******************; |
|---|
| 340 | |
|---|
| 341 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 342 | ;zefield = a_vel_bot |
|---|
| 343 | ;zey = a_h |
|---|
| 344 | ;set_name = 'A_VEL.ps' |
|---|
| 345 | ;set_title = "Scaled vertical velocity variance (m.s!U-1!N)" |
|---|
| 346 | ;set_titlex = 'Scaled vertical velocity variance (m.s!U-1!N)' |
|---|
| 347 | ;set_titley = 'Scaled height (km)' |
|---|
| 348 | ;set_subtitle = '' |
|---|
| 349 | ;set_xrange = [0.,1.2] |
|---|
| 350 | ;set_yrange = [0.,1.2] |
|---|
| 351 | ;set_tickx = 0.1 |
|---|
| 352 | ;set_ticky = 0.1 |
|---|
| 353 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 354 | ;zesim = 1.8 * a_h^(2./3.) * ( 1. - 0.8 * a_h )^2 |
|---|
| 355 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 356 | ;localtimes = localtime[where( (localtime ge 11) and (localtime le 16) )] |
|---|
| 357 | ;;localtimes = localtime |
|---|
| 358 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 359 | ;if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 360 | ;!P.Charsize = 1.2 |
|---|
| 361 | ;altlin=0 |
|---|
| 362 | ;user_lt=localtimes[0] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 363 | ;plot, zefield(*,nntt), zey(*,nntt), xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin, /ISOTROPIC |
|---|
| 364 | ;for ll = 1, n_elements(localtimes)-1 do begin |
|---|
| 365 | ; user_lt=localtimes[ll] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 366 | ; if (nntt ne -1) then oplot, zefield(*,nntt), zey(*,nntt), linestyle=altlin |
|---|
| 367 | ;endfor |
|---|
| 368 | ;loadct, 33 & oplot, zesim, zey, psym=5, color=255 & loadct, 0 |
|---|
| 369 | ;if (saveps eq 'true') then PS_End, /PNG |
|---|
| 370 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 371 | |
|---|
| 372 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 373 | zefield = localtime |
|---|
| 374 | zey = pbl_height |
|---|
| 375 | set_name = 'HEIGHT.ps' |
|---|
| 376 | set_title = 'Boundary layer height (km)' |
|---|
| 377 | set_titlex = 'Local Time (h)' |
|---|
| 378 | set_titley = 'Boundary layer height (km)' |
|---|
| 379 | set_subtitle = '' ;'Criterion is : static stability > 1.5 K.m!U-1!N' |
|---|
| 380 | set_xrange = [8.,17.] |
|---|
| 381 | set_yrange = [0.,7.] |
|---|
| 382 | set_tickx = 1. |
|---|
| 383 | set_ticky = 1. |
|---|
| 384 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 385 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 386 | !P.Charsize = 1.2 |
|---|
| 387 | plot, zefield, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin |
|---|
| 388 | ;oplot, zefield, zey, psym=5 |
|---|
| 389 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 390 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 391 | |
|---|
| 392 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 393 | ;zefield = localtime |
|---|
| 394 | ;zey = w_star_bot |
|---|
| 395 | ;set_name = 'W_STAR.ps' |
|---|
| 396 | ;set_title = "Free convection velocity scale (m.s!U-1!N)" |
|---|
| 397 | ;set_titlex = 'Local Time (h)' |
|---|
| 398 | ;set_titley = 'Free convection velocity scale (m.s!U-1!N)' |
|---|
| 399 | ;set_subtitle = '' |
|---|
| 400 | ;set_xrange = [8.,18.] |
|---|
| 401 | ;set_yrange = [0.,6.] |
|---|
| 402 | ;set_tickx = 1. |
|---|
| 403 | ;set_ticky = 0.5 |
|---|
| 404 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 405 | ;if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 406 | ;!P.Charsize = 1.2 |
|---|
| 407 | ;plot, zefield, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin |
|---|
| 408 | ;;oplot, zefield, zey, psym=5 |
|---|
| 409 | ;if (saveps eq 'true') then PS_End, /PNG |
|---|
| 410 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 411 | |
|---|
| 412 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 413 | ;zefield = localtime |
|---|
| 414 | ;zey = hfbot |
|---|
| 415 | ;zey2 = hfsurf |
|---|
| 416 | ;zey3 = hfbotmax |
|---|
| 417 | ;set_name = 'HFBOT.ps' |
|---|
| 418 | ;set_title = "Heat flux at the top of the radiative layer (K.m.s!U-1!N)" |
|---|
| 419 | ;set_titlex = 'Local Time (h)' |
|---|
| 420 | ;set_titley = 'Heat flux at the top of the radiative layer (K.m.s!U-1!N)' |
|---|
| 421 | ;set_subtitle = '' |
|---|
| 422 | ;set_xrange = [8.,18.] |
|---|
| 423 | ;set_yrange = [0.,2.5] |
|---|
| 424 | ;set_tickx = 1. |
|---|
| 425 | ;set_ticky = 0.5 |
|---|
| 426 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 427 | ;if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 428 | ;!P.Charsize = 1.2 |
|---|
| 429 | ;plot, zefield, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin |
|---|
| 430 | ;oplot, zefield, zey2, linestyle=1 |
|---|
| 431 | ;oplot, zefield, zey3, linestyle=2 |
|---|
| 432 | ;if (saveps eq 'true') then PS_End, /PNG |
|---|
| 433 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 434 | |
|---|
| 435 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 436 | zefield = transpose(tke) |
|---|
| 437 | zex = localtime |
|---|
| 438 | zey = h |
|---|
| 439 | set_name = 'TKE.ps' |
|---|
| 440 | set_title = "Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<u'!U2!N>+<v'!U2!N>+<w'!U2!N>] |
|---|
| 441 | set_titlex = 'Local Time (h)' |
|---|
| 442 | set_titley = 'Altitude above surface (km)' |
|---|
| 443 | set_subtitle = '' ;'Mean over the simulation domain' |
|---|
| 444 | set_xrange = [8.,17.] |
|---|
| 445 | set_yrange = [0.,7.] ;; [0.,200.] |
|---|
| 446 | set_tickx = 1. |
|---|
| 447 | set_ticky = 1. ;; 50. |
|---|
| 448 | minval = 0. |
|---|
| 449 | maxval = 15. |
|---|
| 450 | nlev = maxval-minval |
|---|
| 451 | pal = 22 |
|---|
| 452 | rrr = 'no' |
|---|
| 453 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 454 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 455 | !P.Charsize = 1.2 |
|---|
| 456 | ;; 0. levels |
|---|
| 457 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 458 | ;;; 1. background |
|---|
| 459 | loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 460 | ;; 2. color field |
|---|
| 461 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 462 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 463 | ;; 3. contour field |
|---|
| 464 | loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 465 | ;; 4. choose output |
|---|
| 466 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 467 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 468 | |
|---|
| 469 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 470 | ;zefield = 2.*transpose(ztke) |
|---|
| 471 | ;set_title = "Vertical wind variance (m!U2!N.s!U-2!N)" |
|---|
| 472 | ;minval = -6. |
|---|
| 473 | ;maxval = 12. |
|---|
| 474 | ;pal = 0 |
|---|
| 475 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 476 | zefield = transpose(ztke) |
|---|
| 477 | zex = localtime |
|---|
| 478 | zey = h |
|---|
| 479 | set_name = 'zTKE.ps' |
|---|
| 480 | set_title = "Vertical Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<w'!U2!N>] |
|---|
| 481 | set_titlex = 'Local Time (h)' |
|---|
| 482 | set_titley = 'Altitude above surface (km)' |
|---|
| 483 | set_subtitle = '' ;'Mean over the simulation domain' |
|---|
| 484 | set_xrange = [8.,17.] |
|---|
| 485 | set_yrange = [0.,7.] |
|---|
| 486 | set_tickx = 1. |
|---|
| 487 | set_ticky = 1. |
|---|
| 488 | minval = 0. |
|---|
| 489 | maxval = 8. |
|---|
| 490 | nlev = maxval-minval |
|---|
| 491 | pal = 22 |
|---|
| 492 | rrr = 'no' |
|---|
| 493 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 494 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 495 | !P.Charsize = 1.2 |
|---|
| 496 | ;; 0. levels |
|---|
| 497 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 498 | ;;; 1. background |
|---|
| 499 | loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 500 | ;; 2. color field |
|---|
| 501 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 502 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 503 | ;; 3. contour field |
|---|
| 504 | loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 505 | ;; 4. choose output |
|---|
| 506 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 507 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 508 | |
|---|
| 509 | ;restore, filename='addturb.dat' |
|---|
| 510 | ;modvar = SMOOTH(TEMPORARY(modvar), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 511 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 512 | ;zefield = transpose(modvar) |
|---|
| 513 | ;zex = localtime |
|---|
| 514 | ;zey = h |
|---|
| 515 | ;set_name = 'modvar.ps' |
|---|
| 516 | ;set_title = "Horizontal wind speed variance (m!U2!N.s!U-2!N)" |
|---|
| 517 | ;set_titlex = 'Local Time (h)' |
|---|
| 518 | ;set_titley = 'Altitude above surface (km)' |
|---|
| 519 | ;set_subtitle = '' ;'Mean over the simulation domain' |
|---|
| 520 | ;set_xrange = [7.,17.] |
|---|
| 521 | ;set_yrange = [0.,7.] |
|---|
| 522 | ;set_tickx = 1. |
|---|
| 523 | ;set_ticky = 1. |
|---|
| 524 | ;minval = 0. |
|---|
| 525 | ;maxval = 10. |
|---|
| 526 | ;nlev = maxval-minval |
|---|
| 527 | ;pal = 22 |
|---|
| 528 | ;rrr = 'no' |
|---|
| 529 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 530 | ;if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 531 | ;!P.Charsize = 1.2 |
|---|
| 532 | ;;; 0. levels |
|---|
| 533 | ;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 534 | ;;;; 1. background |
|---|
| 535 | ;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 536 | ;;; 2. color field |
|---|
| 537 | ;loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 538 | ; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 539 | ;;; 3. contour field |
|---|
| 540 | ;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 541 | ;;; 4. choose output |
|---|
| 542 | ;if (saveps eq 'true') then PS_End, /PNG |
|---|
| 543 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 544 | |
|---|
| 545 | |
|---|
| 546 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 547 | zefield = transpose(wt) |
|---|
| 548 | zex = localtime |
|---|
| 549 | zey = h |
|---|
| 550 | set_name = 'HF.ps' |
|---|
| 551 | set_title = "Turbulent Heat Flux (K.m.s!U-1!N)" ;"Vertical Eddy Heat Flux <w'T'>" ;<w'!7h!3'> |
|---|
| 552 | set_titlex = 'Local Time (h)' |
|---|
| 553 | set_titley = 'Altitude above surface (km)' |
|---|
| 554 | set_subtitle = '' ;'Mean over the simulation domain' |
|---|
| 555 | set_xrange = [8.,17.] |
|---|
| 556 | set_yrange = [0.,7.] ;; [0.,200.] |
|---|
| 557 | set_tickx = 1. |
|---|
| 558 | set_ticky = 1. ;; 50. |
|---|
| 559 | minval = -2. |
|---|
| 560 | maxval = 2. |
|---|
| 561 | nlev = floor(maxval-minval)*10 |
|---|
| 562 | pal = 33 |
|---|
| 563 | rrr = 'no' |
|---|
| 564 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 565 | ;minval = -0.8 |
|---|
| 566 | ;maxval = 1.2 |
|---|
| 567 | ;pal = 0 |
|---|
| 568 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 569 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 570 | !P.Charsize = 1.2 |
|---|
| 571 | ;; 0. levels |
|---|
| 572 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 573 | ;;; 1. background |
|---|
| 574 | loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 575 | ;; 2. color field |
|---|
| 576 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 577 | ;;;-------------------------------------------------------------------------------------------------------------------------------- |
|---|
| 578 | ;;; WHITE ZONE - 1. get location of interval in the CT - 2. change the CT to have a white zone |
|---|
| 579 | ulim=0.09 & dlim=-0.09 & w=where(lev le dlim) & n1=w[n_elements(w)-1] & w=where(lev ge ulim) & n2=w[0] & yy=BYTSCL(lev) & nd=yy[n1] & nu=yy[n2]-5 |
|---|
| 580 | nu = nd + (nu-nd)/2 ;; otherwise the interval is too large (because we removed 0) |
|---|
| 581 | TVLCT, r, g, b, /Get & r[nd:nu]=255 & g[nd:nu]=255 & b[nd:nu]=255 & TVLCT, r, g, b |
|---|
| 582 | ;;;-------------------------------------------------------------------------------------------------------------------------------- |
|---|
| 583 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 584 | ;; 3. contour field |
|---|
| 585 | loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 586 | ;; 4. choose output |
|---|
| 587 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 588 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 589 | |
|---|
| 590 | ;;restore, filename='tpot_profB' |
|---|
| 591 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 592 | ;zefield = t |
|---|
| 593 | ;zey = ht |
|---|
| 594 | ;set_name = 'T.ps' |
|---|
| 595 | ;set_title = "Potential Temperature (K)" |
|---|
| 596 | ;set_titlex = 'Potential Temperature (K)' |
|---|
| 597 | ;set_titley = 'Altitude above surface (km)' |
|---|
| 598 | ;set_subtitle = '' ;'Mean over the simulation domain' |
|---|
| 599 | ;set_xrange = [190.,230.] ;[min(t),max(t)] |
|---|
| 600 | ;set_yrange = [0.,7.] |
|---|
| 601 | ;set_tickx = 5. |
|---|
| 602 | ;set_ticky = 1. |
|---|
| 603 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 604 | ;;localtimes = [9,10,11,12,13,14,15,16,17,18] |
|---|
| 605 | ;localtimes = [7,9,11,13,15,17] |
|---|
| 606 | ;;localtimes = [17.0,17.1,17.2,17.3,17.4,17.5,17.6,17.7,17.8,17.9,18.0] |
|---|
| 607 | ;;localtimes = [15.0] |
|---|
| 608 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 609 | ;if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 610 | ;!P.Charsize = 1.2 |
|---|
| 611 | ;;altlin=0 & loadct, 0 |
|---|
| 612 | ;altlin=1 & loadct, 0 |
|---|
| 613 | ;user_lt=localtimes[0] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 614 | ;plot, zefield(*,nntt), zey(*,nntt), xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin |
|---|
| 615 | ;;plot, zefield(*,nntt), les_column, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin |
|---|
| 616 | ;for ll = 1, n_elements(localtimes)-1 do begin |
|---|
| 617 | ; CASE altlin OF |
|---|
| 618 | ; 0: altlin=1 |
|---|
| 619 | ; 1: altlin=0 |
|---|
| 620 | ; ENDCASE |
|---|
| 621 | ; user_lt=localtimes[ll] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 622 | ; if (nntt ne -1) then oplot, zefield(*,nntt), zey(*,nntt), linestyle=altlin |
|---|
| 623 | ;; if (nntt ne -1) then oplot, zefield(*,nntt), les_column, linestyle=altlin |
|---|
| 624 | ;endfor |
|---|
| 625 | ;;oplot, ro_tpot, ro_column, psym=7 |
|---|
| 626 | ;;xyouts, 192.0, 0.25, '09:00' |
|---|
| 627 | ;;xyouts, 205.5, 0.25, '11:00' |
|---|
| 628 | ;;xyouts, 214.5, 0.25, '13:00' |
|---|
| 629 | ;;xyouts, 219.5, 0.25, '17:00' |
|---|
| 630 | ;;xyouts, 220.5, 2.30, '17:00 [RO]' |
|---|
| 631 | ;if (saveps eq 'true') then PS_End, /PNG |
|---|
| 632 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 633 | |
|---|
| 634 | |
|---|
| 635 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 636 | ;zefield = stst |
|---|
| 637 | ;zey = ht |
|---|
| 638 | ;set_name = 'STST.ps' |
|---|
| 639 | ;set_title = "Static stability (K.m!U-1!N)" |
|---|
| 640 | ;set_titlex = 'Static stability (K.m!U-1!N)' |
|---|
| 641 | ;set_titley = 'Altitude above surface (km)' |
|---|
| 642 | ;set_subtitle = 'Mean over the simulation domain' |
|---|
| 643 | ;set_xrange = [-5.,5.] |
|---|
| 644 | ;set_yrange = [0.,7.] |
|---|
| 645 | ;set_tickx = 1. |
|---|
| 646 | ;set_ticky = 1. |
|---|
| 647 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 648 | ;localtimes = [9,11,13,15,17] |
|---|
| 649 | ;localtimes = [12,13,14,15,16] |
|---|
| 650 | ;;localtimes = [17.0,17.1,17.2,17.3,17.4,17.5,17.6,17.7,17.8,17.9,18.0] |
|---|
| 651 | ;;localtimes = [15.0] |
|---|
| 652 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 653 | ;if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 654 | ;!P.Charsize = 1.2 |
|---|
| 655 | ;altlin=0 & loadct, 0 |
|---|
| 656 | ;user_lt=localtimes[0] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 657 | ;plot, zefield(*,nntt), zey(*,nntt), xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin |
|---|
| 658 | ;oplot, zefield(*,nntt), zey(*,nntt), psym=5 |
|---|
| 659 | ;for ll = 1, n_elements(localtimes)-1 do begin |
|---|
| 660 | ; CASE altlin OF |
|---|
| 661 | ; 0: altlin=1 |
|---|
| 662 | ; 1: altlin=0 |
|---|
| 663 | ; ENDCASE |
|---|
| 664 | ; user_lt=localtimes[ll] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 665 | ; if (nntt ne -1) then oplot, zefield(*,nntt), zey(*,nntt), linestyle=altlin |
|---|
| 666 | ; if (nntt ne -1) then oplot, zefield(*,nntt), zey(*,nntt), psym=5 |
|---|
| 667 | ;endfor |
|---|
| 668 | ;oplot, 0.*zefield(*,nntt) + 1.5, zey(*,nntt), linestyle=2 |
|---|
| 669 | ;if (saveps eq 'true') then PS_End, /PNG |
|---|
| 670 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 671 | |
|---|
| 672 | |
|---|
| 673 | |
|---|
| 674 | if (n_elements(xtke) eq 0) then stop |
|---|
| 675 | |
|---|
| 676 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 677 | ;;zefield = 2.*transpose(xtke) |
|---|
| 678 | ;;minval = -4 |
|---|
| 679 | ;;maxval = 8. ;12. |
|---|
| 680 | ;;pal = 0 |
|---|
| 681 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 682 | ;zex = localtime |
|---|
| 683 | ;zey = h |
|---|
| 684 | ;set_name = 'xTKE.ps' |
|---|
| 685 | ;set_title = "Horizontal Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<u'!U2!N>] |
|---|
| 686 | ;set_titlex = 'Local Time (h)' |
|---|
| 687 | ;set_titley = 'Altitude above surface (km)' |
|---|
| 688 | ;set_subtitle = '';'Mean over the simulation domain' |
|---|
| 689 | ;set_xrange = [8.,17.] |
|---|
| 690 | ;set_yrange = [0.,7.] |
|---|
| 691 | ;set_tickx = 1. |
|---|
| 692 | ;set_ticky = 1. |
|---|
| 693 | ;minval = 0. |
|---|
| 694 | ;maxval = 8. |
|---|
| 695 | ;nlev = maxval-minval |
|---|
| 696 | ;pal = 22 |
|---|
| 697 | ;rrr = 'no' |
|---|
| 698 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 699 | ;if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 700 | ;!P.Charsize = 1.2 |
|---|
| 701 | ;;; 0. levels |
|---|
| 702 | ;lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 703 | ;;;; 1. background |
|---|
| 704 | ;loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 705 | ;;; 2. color field |
|---|
| 706 | ;loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 707 | ; contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 708 | ;;; 3. contour field |
|---|
| 709 | ;loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 710 | ;;; 4. choose output |
|---|
| 711 | ;if (saveps eq 'true') then PS_End, /PNG |
|---|
| 712 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 713 | |
|---|
| 714 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 715 | zefield = transpose(ytke) + transpose(xtke) |
|---|
| 716 | zex = localtime |
|---|
| 717 | zey = h |
|---|
| 718 | set_name = 'hTKE.ps' |
|---|
| 719 | set_title = "Horizontal Turbulent Kinetic Energy (m!U2!N.s!U-2!N)" ;0.5[<v'!U2!N>] |
|---|
| 720 | set_titlex = 'Local Time (h)' |
|---|
| 721 | set_titley = 'Altitude above surface (km)' |
|---|
| 722 | set_subtitle = '' ;'Mean over the simulation domain' |
|---|
| 723 | set_xrange = [8.,17.] |
|---|
| 724 | set_yrange = [0.,1.] ;[0.,7.] ;; [0.,200.] |
|---|
| 725 | set_tickx = 1. |
|---|
| 726 | set_ticky = 0.1 ;1. ;50. |
|---|
| 727 | minval = 0. |
|---|
| 728 | maxval = 12. ;10. |
|---|
| 729 | nlev = maxval-minval |
|---|
| 730 | pal = 22 |
|---|
| 731 | rrr = 'no' |
|---|
| 732 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 733 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 734 | !P.Charsize = 1.2 |
|---|
| 735 | ;; 0. levels |
|---|
| 736 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 737 | ;;; 1. background |
|---|
| 738 | loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 739 | ;; 2. color field |
|---|
| 740 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 741 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 742 | ;; 3. contour field |
|---|
| 743 | loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 744 | ;; 4. choose output |
|---|
| 745 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 746 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 747 | |
|---|
| 748 | if (n_elements(wmax) eq 0) then stop |
|---|
| 749 | |
|---|
| 750 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 751 | zefield = transpose(wmax) |
|---|
| 752 | zex = localtime |
|---|
| 753 | zey = h |
|---|
| 754 | set_name = 'Wmax.ps' |
|---|
| 755 | set_title = "Maximum updraft speed (m.s!U-1!N)" |
|---|
| 756 | set_titlex = 'Local Time (h)' |
|---|
| 757 | set_titley = 'Altitude above surface (km)' |
|---|
| 758 | set_subtitle = '' |
|---|
| 759 | set_xrange = [8.,18.] |
|---|
| 760 | set_yrange = [0.,7.] |
|---|
| 761 | set_tickx = 1. |
|---|
| 762 | set_ticky = 1. |
|---|
| 763 | minval = 0. |
|---|
| 764 | maxval = 16. ;15.;14.;13. |
|---|
| 765 | nlev = maxval-minval |
|---|
| 766 | pal = 22 |
|---|
| 767 | rrr = 'no' |
|---|
| 768 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 769 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 770 | !P.Charsize = 1.2 |
|---|
| 771 | ;; 0. levels |
|---|
| 772 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 773 | ;;; 1. background |
|---|
| 774 | loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 775 | ;; 2. color field |
|---|
| 776 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 777 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 778 | ;; 3. contour field |
|---|
| 779 | loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 780 | ;; 4. choose output |
|---|
| 781 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 782 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 783 | |
|---|
| 784 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 785 | zefield = abs(transpose(wmin)) |
|---|
| 786 | zex = localtime |
|---|
| 787 | zey = h |
|---|
| 788 | set_name = 'Wmin.ps' |
|---|
| 789 | set_title = "Maximum downdraft speed (m.s!U-1!N)" |
|---|
| 790 | set_titlex = 'Local Time (h)' |
|---|
| 791 | set_titley = 'Altitude above surface (km)' |
|---|
| 792 | set_subtitle = '' |
|---|
| 793 | set_xrange = [8.,18.] |
|---|
| 794 | set_yrange = [0.,7.] |
|---|
| 795 | set_tickx = 1. |
|---|
| 796 | set_ticky = 1. |
|---|
| 797 | minval = 0. |
|---|
| 798 | maxval = 9.;12. |
|---|
| 799 | nlev = maxval-minval |
|---|
| 800 | pal = 22 |
|---|
| 801 | rrr = 'no' ;'yes' |
|---|
| 802 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 803 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 804 | !P.Charsize = 1.2 |
|---|
| 805 | ;; 0. levels |
|---|
| 806 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 807 | ;;; 1. background |
|---|
| 808 | loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 809 | ;; 2. color field |
|---|
| 810 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 811 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 812 | ;; 3. contour field |
|---|
| 813 | loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 814 | ;; 4. choose output |
|---|
| 815 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 816 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 817 | |
|---|
| 818 | |
|---|
| 819 | restore, filename='addturb.dat' |
|---|
| 820 | velmax = SMOOTH(TEMPORARY(velmax), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 821 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 822 | zefield = transpose(velmax) |
|---|
| 823 | zex = localtime |
|---|
| 824 | zey = h |
|---|
| 825 | set_name = 'velmax.ps' |
|---|
| 826 | set_title = "Maximum horizontal wind speed (m.s!U-1!N)" |
|---|
| 827 | set_titlex = 'Local Time (h)' |
|---|
| 828 | set_titley = 'Altitude above surface (km)' |
|---|
| 829 | set_subtitle = '' ;'Mean over the simulation domain' |
|---|
| 830 | set_xrange = [8.,17.] |
|---|
| 831 | set_yrange = [0.,4.] |
|---|
| 832 | set_tickx = 1. |
|---|
| 833 | set_ticky = 1. |
|---|
| 834 | minval = 0. |
|---|
| 835 | maxval = 30. ;12. |
|---|
| 836 | nlev = maxval-minval |
|---|
| 837 | pal = 22 |
|---|
| 838 | rrr = 'no' |
|---|
| 839 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 840 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 841 | !P.Charsize = 1.2 |
|---|
| 842 | ;; 0. levels |
|---|
| 843 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 844 | ;;; 1. background |
|---|
| 845 | loadct, 0 & contour, /NODATA, zefield, zex, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0 |
|---|
| 846 | ;;; 2. color field |
|---|
| 847 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 848 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 849 | ;; 3. contour field |
|---|
| 850 | loadct, 0 & contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /noerase, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, C_LINESTYLE = (lev LT 0.0) |
|---|
| 851 | ;; 4. choose output |
|---|
| 852 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 853 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 854 | |
|---|
| 855 | |
|---|
| 856 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 857 | ;zefield = localtime |
|---|
| 858 | ;zey = depressions |
|---|
| 859 | ;set_name = 'DEPR.ps' |
|---|
| 860 | ;set_title = "Percentage of depressions in the domain" |
|---|
| 861 | ;set_titlex = 'Local Time (h)' |
|---|
| 862 | ;set_titley = 'Percentage' |
|---|
| 863 | ;set_subtitle = '' |
|---|
| 864 | ;set_xrange = [8.,18.] |
|---|
| 865 | ;set_yrange = [0.,3.] |
|---|
| 866 | ;set_tickx = 1. |
|---|
| 867 | ;set_ticky = 1. |
|---|
| 868 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 869 | ;if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 870 | ;!P.Charsize = 1.2 |
|---|
| 871 | ;plot, zefield, zey, xtitle=set_titlex, xrange=set_xrange, xtickinterval=set_tickx, ytitle=set_titley, yrange=set_yrange, ytickinterval=set_ticky, title=set_title, subtitle=set_subtitle, color=0, linestyle=altlin |
|---|
| 872 | ;oplot, zefield, zey, psym=5 |
|---|
| 873 | ; oplot, localtime, abs(psmin), linestyle=1 |
|---|
| 874 | ;if (saveps eq 'true') then PS_End, /PNG |
|---|
| 875 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 876 | |
|---|
| 877 | stop |
|---|
| 878 | |
|---|
| 879 | |
|---|
| 880 | |
|---|
| 881 | |
|---|
| 882 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 883 | goto, no_staticstab |
|---|
| 884 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 885 | if (saveps eq 'true') then begin |
|---|
| 886 | device, /close |
|---|
| 887 | set_plot, 'ps' & device, filename='plot/staticstab.ps' |
|---|
| 888 | endif |
|---|
| 889 | |
|---|
| 890 | user_lt=17. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) |
|---|
| 891 | |
|---|
| 892 | ;;; recompute height @ given local time |
|---|
| 893 | caca=heightp+hgtu/1000. |
|---|
| 894 | height=reform(ph(*,0,*,*)) |
|---|
| 895 | height=total(height,1)/n_elements(height(*,0,0)) |
|---|
| 896 | height=reform(height) |
|---|
| 897 | height=reform(height(*,yeah(0))) |
|---|
| 898 | height=height/1000./3.72 |
|---|
| 899 | heightp=height(0:n_elements(height(*))-2) |
|---|
| 900 | print, 'new minus old', heightp - caca |
|---|
| 901 | ;;; recompute height @ given local time |
|---|
| 902 | |
|---|
| 903 | press=reform(p(*,0,*,*)) |
|---|
| 904 | press=total(press,1)/n_elements(press(*,0,0)) |
|---|
| 905 | press=reform(press) |
|---|
| 906 | press=reform(press(*,yeah(0))) |
|---|
| 907 | press=press(0:n_elements(press(*))-2) |
|---|
| 908 | |
|---|
| 909 | staticstab = DERIV( heightp , reform(what_I_plot5(*,yeah(0))) ) + 1000.*3.72/844.6 ;; 4.9 dans Hinson |
|---|
| 910 | staticstab = staticstab(0:n_elements(staticstab)-5) |
|---|
| 911 | heightp = heightp(0:n_elements(heightp)-5) |
|---|
| 912 | |
|---|
| 913 | plot, $ |
|---|
| 914 | staticstab, $ |
|---|
| 915 | heightp,$ |
|---|
| 916 | xtitle='Static stability (K/km)',$ |
|---|
| 917 | xrange=[-1.,5.], $ |
|---|
| 918 | xtickinterval=0.5, $ |
|---|
| 919 | ytitle='Altitude above surface (km)', $ |
|---|
| 920 | yrange=[min(heightp),max(heightp)], $ |
|---|
| 921 | ytickinterval=1., $ |
|---|
| 922 | title="LMD LES Static stability @ LT 17h (K/km)", $ |
|---|
| 923 | subtitle='zonal average at lat. '+latwrite |
|---|
| 924 | oplot, $ |
|---|
| 925 | staticstab, $ |
|---|
| 926 | heightp,$ |
|---|
| 927 | psym=5 |
|---|
| 928 | oplot, $ |
|---|
| 929 | findgen(n_elements(heightp))*0. + 1.,$ |
|---|
| 930 | heightp,$ |
|---|
| 931 | linestyle=2 |
|---|
| 932 | oplot, $ |
|---|
| 933 | findgen(n_elements(heightp))*0. + 2.,$ |
|---|
| 934 | heightp,$ |
|---|
| 935 | linestyle=2 |
|---|
| 936 | |
|---|
| 937 | ;;;;;;;;;; |
|---|
| 938 | w = where((staticstab gt 1.5) and ((heightp- hgtu/1000.) gt 1.)) |
|---|
| 939 | t_top_plus = what_I_plot5(w(0),yeah(0)) |
|---|
| 940 | z_top_plus = heightp(w(0)) |
|---|
| 941 | p_top_plus = press(w(0)) |
|---|
| 942 | t_top_moins = what_I_plot5(w(0)-1,yeah(0)) |
|---|
| 943 | z_top_moins = heightp(w(0)-1) |
|---|
| 944 | p_top_moins = press(w(0)-1) |
|---|
| 945 | pbl_depth = (z_top_plus*alog(p_top_plus) + z_top_moins*alog(p_top_moins))/(alog(p_top_plus) + alog(p_top_moins)) - hgtu/1000. |
|---|
| 946 | xyouts, 3., 1.5 + (max(heightp) + min(heightp)) / 3., 'Ls = '+string(lsu,'(F5.1)')+'!Uo!N', CHARSIZE=1 |
|---|
| 947 | xyouts, 3., 1. + (max(heightp) + min(heightp)) / 3., 'Lat = '+string(latu,'(F5.1)')+'!Uo!N', CHARSIZE=1 |
|---|
| 948 | xyouts, 3., 0.5 + (max(heightp) + min(heightp)) / 3., 'LonE = '+string(lonu,'(F6.1)')+'!Uo!N', CHARSIZE=1 |
|---|
| 949 | xyouts, 3., (max(heightp) + min(heightp)) / 3., 'T!Dt!N = '+string(t_top_plus,'(I0)')+'/'+string(t_top_moins,'(I0)')+' K ', CHARSIZE=1 |
|---|
| 950 | xyouts, 3., -0.5 + (max(heightp) + min(heightp)) / 3., 'p!Dt!N = '+string(p_top_plus,'(I0)')+'/'+string(p_top_moins,'(I0)')+' Pa', CHARSIZE=1 |
|---|
| 951 | xyouts, 3., -1. + (max(heightp) + min(heightp)) / 3., 'z!Dt!N = '+string(z_top_plus,'(F4.1)')+'/'+string(z_top_moins,'(F4.1)')+' km', CHARSIZE=1 |
|---|
| 952 | xyouts, 3., -1.5 + (max(heightp) + min(heightp)) / 3., 'z!Ds!N = '+string(hgtu/1000.,'(F4.1)')+' km', CHARSIZE=1 |
|---|
| 953 | xyouts, 3., -2. + (max(heightp) + min(heightp)) / 3., 'D = '+string(pbl_depth,'(F4.1)')+' km', CHARSIZE=1 |
|---|
| 954 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 955 | no_staticstab: |
|---|
| 956 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 957 | |
|---|
| 958 | |
|---|
| 959 | |
|---|
| 960 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 961 | if (saveps eq 'true') then begin |
|---|
| 962 | device, /close |
|---|
| 963 | endif |
|---|
| 964 | stop |
|---|
| 965 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 966 | |
|---|
| 967 | |
|---|
| 968 | user_lt=13. |
|---|
| 969 | yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) |
|---|
| 970 | |
|---|
| 971 | user_h=0.1 |
|---|
| 972 | walt=where(abs(heightp-user_h) eq (min(abs(heightp-user_h)))) |
|---|
| 973 | |
|---|
| 974 | ;mapfield=reform(t[*,*,walt(0),w(0)]) |
|---|
| 975 | ;help, mapfield |
|---|
| 976 | ;contour, mapfield |
|---|
| 977 | |
|---|
| 978 | section=reform(w[*,0,*,yeah(0)]) |
|---|
| 979 | |
|---|
| 980 | section=reform(w[*,0,*,160]) |
|---|
| 981 | |
|---|
| 982 | lev=[-12.,-8.,-4.,4.,8.,12.] |
|---|
| 983 | lev=[-5.,-4.,-3.,-2.,-1.,1.,2.,3.,4.,5.] |
|---|
| 984 | contour, $ |
|---|
| 985 | section, $ |
|---|
| 986 | (lon(*,0)-lon(0,0))*59., $ |
|---|
| 987 | heightp, $ |
|---|
| 988 | levels=lev , $ |
|---|
| 989 | C_LINESTYLE = (lev LT 0.0) |
|---|
| 990 | |
|---|
| 991 | device, /close |
|---|
| 992 | |
|---|
| 993 | end |
|---|
| 994 | |
|---|
