| 1 | pro getturb, saveps=saveps |
|---|
| 2 | |
|---|
| 3 | |
|---|
| 4 | |
|---|
| 5 | ;---------------------------------- |
|---|
| 6 | ; USE: getturb |
|---|
| 7 | ; getturb, saveps='false' |
|---|
| 8 | ;---------------------------------- |
|---|
| 9 | |
|---|
| 10 | |
|---|
| 11 | history_interval_s = 100. |
|---|
| 12 | smoothampl=3700/history_interval_s |
|---|
| 13 | smoothampl=0. |
|---|
| 14 | |
|---|
| 15 | ; |
|---|
| 16 | ; constantes |
|---|
| 17 | ; |
|---|
| 18 | p0=610. & t0=220. & r_cp=1/4.4 & grav=3.72 & R=192. |
|---|
| 19 | |
|---|
| 20 | ; |
|---|
| 21 | ; graphics definition |
|---|
| 22 | ; |
|---|
| 23 | if (n_elements(saveps) eq 0) then saveps='true' |
|---|
| 24 | if (saveps eq 'false') then begin |
|---|
| 25 | ;!p.multi=[0,3,2] |
|---|
| 26 | !P.CHARSIZE=2. |
|---|
| 27 | WINDOW, /PIXMAP & WDELETE & DEVICE,BYPASS_TRANSLATION=0,DECOMPOSED=0,RETAIN=2 |
|---|
| 28 | endif else begin |
|---|
| 29 | PREF_SET, 'IDL_PATH', '/home/spiga/Save/SOURCES/IDL/fsc_psconfig:<IDL_DEFAULT>', /COMMIT |
|---|
| 30 | endelse |
|---|
| 31 | |
|---|
| 32 | ; |
|---|
| 33 | ; retrieve fields |
|---|
| 34 | ; |
|---|
| 35 | openr,unit,'getturb.dat',/get_lun,error=err |
|---|
| 36 | IF (err ne 0) THEN BEGIN |
|---|
| 37 | |
|---|
| 38 | ; |
|---|
| 39 | ; input files |
|---|
| 40 | ; |
|---|
| 41 | OPENR, 22, 'input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22 |
|---|
| 42 | OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23 |
|---|
| 43 | |
|---|
| 44 | ; |
|---|
| 45 | ; get fields |
|---|
| 46 | ; |
|---|
| 47 | domain='d01' & filesWRF = FindFile('wrfout_'+domain+'_????-??-??_??:??:??') & nf=n_elements(filesWRF) |
|---|
| 48 | |
|---|
| 49 | ; |
|---|
| 50 | ; get dimensions |
|---|
| 51 | ; |
|---|
| 52 | id=ncdf_open(filesWRF(0)) |
|---|
| 53 | NCDF_DIMINQ, id, NCDF_DIMID(id, 'west_east' ), toto, nx & NCDF_DIMINQ, id, NCDF_DIMID(id, 'south_north' ), toto, ny |
|---|
| 54 | NCDF_DIMINQ, id, NCDF_DIMID(id, 'bottom_top' ), toto, nz & NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), toto, nt |
|---|
| 55 | NCDF_CLOSE, id |
|---|
| 56 | |
|---|
| 57 | ; |
|---|
| 58 | ; prepare loop |
|---|
| 59 | ; |
|---|
| 60 | nloop1 = nf-1 & nloop2 = nt & yeye = 0 |
|---|
| 61 | localtime = lctu + history_interval_s*findgen(nloop1*nloop2)/3700. |
|---|
| 62 | wt = fltarr(nz,nloop1*nloop2) & tke = fltarr(nz,nloop1*nloop2) & ztke = fltarr(nz,nloop1*nloop2) & t = fltarr(nz,nloop1*nloop2) |
|---|
| 63 | p = fltarr(nz) & ph = fltarr(nz) & pht = fltarr(nz,nloop1*nloop2) & pt = fltarr(nz,nloop1*nloop2) & stst = fltarr(nz,nloop1*nloop2) |
|---|
| 64 | |
|---|
| 65 | ; |
|---|
| 66 | ; loop loop |
|---|
| 67 | ; |
|---|
| 68 | for loop = 0, nloop1-1 do begin |
|---|
| 69 | timetime = SYSTIME(1) |
|---|
| 70 | for loop2 = 0, nloop2-1 do begin |
|---|
| 71 | |
|---|
| 72 | ; t' = t - <t> |
|---|
| 73 | ; ------------ |
|---|
| 74 | yeyeye = 1. |
|---|
| 75 | tprime = getget(filesWRF(loop), 'T', anomaly=yeyeye, count=[0,0,0,1], offset=[0,0,0,loop2]) ;; t' = t - <t> |
|---|
| 76 | t(*,yeye) = t0 + TEMPORARY(yeyeye) |
|---|
| 77 | ; w' = w |
|---|
| 78 | ; ------ |
|---|
| 79 | wprime = getget(filesWRF(loop), 'W', count=[0,0,0,1], offset=[0,0,0,loop2]) |
|---|
| 80 | ; tke = 0.5 ( <u'^2> + <v'^2> + <w'^2> ) ; u' = u ; v' = v |
|---|
| 81 | ; -------------------------------------------------------- |
|---|
| 82 | ztke(*,yeye) = 0.5 * TOTAL(TOTAL(wprime^2,1),1) / float(nx) / float(ny) |
|---|
| 83 | tke(*,yeye) = ztke(*,yeye) + $ |
|---|
| 84 | 0.5 * ( $ |
|---|
| 85 | TOTAL(TOTAL(getget(filesWRF(loop), 'U', count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) + $ |
|---|
| 86 | TOTAL(TOTAL(getget(filesWRF(loop), 'V', count=[0,0,0,1], offset=[0,0,0,loop2])^2,1),1) $ |
|---|
| 87 | ) / float(nx) / float(ny) |
|---|
| 88 | ; <w't'> |
|---|
| 89 | ; ------ |
|---|
| 90 | wt(*,yeye) = TOTAL(TOTAL(TEMPORARY(tprime) * TEMPORARY(wprime),1),1) / float(nx) / float(ny) |
|---|
| 91 | ; p & ph |
|---|
| 92 | ; ------ |
|---|
| 93 | if (loop + loop2 eq 0) then nloopbeware = nloop2-1 else nloopbeware = nloop2 ; 1ere valeur vaut 0 |
|---|
| 94 | 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 |
|---|
| 95 | pt(*,yeye) = TOTAL(TOTAL(getget(filesWRF(loop), 'PTOT' , count=[0,0,0,1], offset=[0,0,0,loop2]),1),1) / float(nx) / float(ny) |
|---|
| 96 | ph = TEMPORARY(ph) + pht(*,yeye) / nloopbeware / nloop1 |
|---|
| 97 | p = TEMPORARY(p ) + pt(*,yeye) / nloop2 / nloop1 |
|---|
| 98 | ; static stability |
|---|
| 99 | ; ---------------- |
|---|
| 100 | 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 |
|---|
| 101 | ; loop count |
|---|
| 102 | ; --------- |
|---|
| 103 | yeye = TEMPORARY(yeye) + 1 |
|---|
| 104 | |
|---|
| 105 | endfor |
|---|
| 106 | print, 'file '+string(loop+1,'(I0)'), SYSTIME(1) - timetime, ' s' |
|---|
| 107 | endfor |
|---|
| 108 | h = TEMPORARY(ph) - hgtu/1000. ;; altitude above ground |
|---|
| 109 | ht = TEMPORARY(pht) - hgtu/1000. |
|---|
| 110 | ; |
|---|
| 111 | ; save |
|---|
| 112 | ; |
|---|
| 113 | save, wt, tke, ztke, h, ht, t, p, pt, stst, localtime, filename='getturb.dat' |
|---|
| 114 | |
|---|
| 115 | ENDIF ELSE BEGIN |
|---|
| 116 | |
|---|
| 117 | print, 'OK, file is here' |
|---|
| 118 | restore, filename='getturb.dat' |
|---|
| 119 | |
|---|
| 120 | ENDELSE |
|---|
| 121 | |
|---|
| 122 | ; |
|---|
| 123 | ; smooth smooth |
|---|
| 124 | ; |
|---|
| 125 | wt = SMOOTH(TEMPORARY(wt), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 126 | tke = SMOOTH(TEMPORARY(tke), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 127 | ztke = SMOOTH(TEMPORARY(ztke), [0,smoothampl], /EDGE_TRUNCATE) |
|---|
| 128 | |
|---|
| 129 | |
|---|
| 130 | ;*******************; |
|---|
| 131 | ;*******************; |
|---|
| 132 | ; PLOTS PLOTS PLOTS ; |
|---|
| 133 | ;*******************; |
|---|
| 134 | ;*******************; |
|---|
| 135 | |
|---|
| 136 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 137 | zefield = transpose(tke) |
|---|
| 138 | zex = localtime |
|---|
| 139 | zey = h |
|---|
| 140 | set_name = 'TKE.ps' |
|---|
| 141 | set_title = "Turbulent Kinetic Energy 0.5[<u'!U2!N>+<v'!U2!N>+<w'!U2!N>] (m!U2!N.s!U-2!N)" |
|---|
| 142 | set_titlex = 'Local Time (h)' |
|---|
| 143 | set_titley = 'Altitude above surface (km)' |
|---|
| 144 | set_subtitle = 'Mean over the simulation domain' |
|---|
| 145 | set_xrange = [8.,18.] |
|---|
| 146 | set_yrange = [0.,10.] |
|---|
| 147 | set_tickx = 1. |
|---|
| 148 | set_ticky = 1. |
|---|
| 149 | minval = 0. |
|---|
| 150 | maxval = 20. ;15. |
|---|
| 151 | nlev = maxval-minval |
|---|
| 152 | pal = 22 |
|---|
| 153 | rrr = 'no' |
|---|
| 154 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 155 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 156 | !P.Charsize = 1.2 |
|---|
| 157 | ;; 0. levels |
|---|
| 158 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 159 | ;;; 1. background |
|---|
| 160 | 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 |
|---|
| 161 | ;; 2. color field |
|---|
| 162 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 163 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 164 | ;; 3. contour field |
|---|
| 165 | 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) |
|---|
| 166 | ;; 4. choose output |
|---|
| 167 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 168 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 169 | |
|---|
| 170 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 171 | zefield = transpose(ztke) |
|---|
| 172 | ;zefield = 100.*transpose(ztke)/transpose(tke) & zefield[where(transpose(tke) lt 1.)]=0. |
|---|
| 173 | zex = localtime |
|---|
| 174 | zey = h |
|---|
| 175 | set_name = 'zTKE.ps' |
|---|
| 176 | set_title = "Vertical Turbulent Kinetic Energy 0.5[<w'!U2!N>] (m!U2!N.s!U-2!N)" |
|---|
| 177 | set_titlex = 'Local Time (h)' |
|---|
| 178 | set_titley = 'Altitude above surface (km)' |
|---|
| 179 | set_subtitle = 'Mean over the simulation domain' |
|---|
| 180 | set_xrange = [8.,18.] |
|---|
| 181 | set_yrange = [0.,10.] |
|---|
| 182 | set_tickx = 1. |
|---|
| 183 | set_ticky = 1. |
|---|
| 184 | minval = 0. |
|---|
| 185 | maxval = 14. ;10. |
|---|
| 186 | nlev = maxval-minval |
|---|
| 187 | pal = 22 |
|---|
| 188 | rrr = 'no' |
|---|
| 189 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 190 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 191 | !P.Charsize = 1.2 |
|---|
| 192 | ;; 0. levels |
|---|
| 193 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 194 | ;;; 1. background |
|---|
| 195 | 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 |
|---|
| 196 | ;; 2. color field |
|---|
| 197 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 198 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 199 | ;; 3. contour field |
|---|
| 200 | 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) |
|---|
| 201 | ;; 4. choose output |
|---|
| 202 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 203 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 204 | |
|---|
| 205 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 206 | zefield = transpose(wt) |
|---|
| 207 | zex = localtime |
|---|
| 208 | zey = h |
|---|
| 209 | set_name = 'HF.ps' |
|---|
| 210 | set_title = "Vertical Eddy Heat Flux <w'!7h!3'> (K.m.s!U-1!N)" |
|---|
| 211 | set_titlex = 'Local Time (h)' |
|---|
| 212 | set_titley = 'Altitude above surface (km)' |
|---|
| 213 | set_subtitle = 'Mean over the simulation domain' |
|---|
| 214 | set_xrange = [8.,18.] |
|---|
| 215 | set_yrange = [0.,10.] |
|---|
| 216 | set_tickx = 1. |
|---|
| 217 | set_ticky = 1. |
|---|
| 218 | minval = -2. |
|---|
| 219 | maxval = 2. |
|---|
| 220 | nlev = floor(maxval-minval)*10 |
|---|
| 221 | pal = 33 ;4 |
|---|
| 222 | rrr = 'no' |
|---|
| 223 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 224 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 225 | !P.Charsize = 1.2 |
|---|
| 226 | ;; 0. levels |
|---|
| 227 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 228 | ;;; 1. background |
|---|
| 229 | 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 |
|---|
| 230 | ;; 2. color field |
|---|
| 231 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 232 | ;;;-------------------------------------------------------------------------------------------------------------------------------- |
|---|
| 233 | ;;; WHITE ZONE - 1. get location of interval in the CT - 2. change the CT to have a white zone |
|---|
| 234 | 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 |
|---|
| 235 | nu = nd + (nu-nd)/2 ;; otherwise the interval is too large (because we removed 0) |
|---|
| 236 | TVLCT, r, g, b, /Get & r[nd:nu]=255 & g[nd:nu]=255 & b[nd:nu]=255 & TVLCT, r, g, b |
|---|
| 237 | ;;;-------------------------------------------------------------------------------------------------------------------------------- |
|---|
| 238 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 239 | ;; 3. contour field |
|---|
| 240 | 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) |
|---|
| 241 | ;; 4. choose output |
|---|
| 242 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 243 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 244 | |
|---|
| 245 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 246 | zefield = t |
|---|
| 247 | zey = ht |
|---|
| 248 | set_name = 'T.ps' |
|---|
| 249 | set_title = "Potential Temperature (K)" |
|---|
| 250 | set_titlex = 'Potential Temperature (K)' |
|---|
| 251 | set_titley = 'Altitude above surface (km)' |
|---|
| 252 | set_subtitle = 'Mean over the simulation domain' |
|---|
| 253 | set_xrange = [min(t),max(t)] |
|---|
| 254 | set_yrange = [0.,10.] |
|---|
| 255 | set_tickx = 5. |
|---|
| 256 | set_ticky = 1. |
|---|
| 257 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 258 | localtimes = [9,10,11,12,13,14,15,16,17,18] |
|---|
| 259 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 260 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 261 | !P.Charsize = 1.2 |
|---|
| 262 | altlin=0 & loadct, 0 |
|---|
| 263 | user_lt=localtimes[0] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 264 | 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 |
|---|
| 265 | for ll = 1, n_elements(localtimes)-1 do begin |
|---|
| 266 | CASE altlin OF |
|---|
| 267 | 0: altlin=1 |
|---|
| 268 | 1: altlin=0 |
|---|
| 269 | ENDCASE |
|---|
| 270 | user_lt=localtimes[ll] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 271 | if (nntt ne -1) then oplot, zefield(*,nntt), zey(*,nntt), linestyle=altlin |
|---|
| 272 | endfor |
|---|
| 273 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 274 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 275 | |
|---|
| 276 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 277 | zefield = stst |
|---|
| 278 | zey = ht |
|---|
| 279 | set_name = 'STST.ps' |
|---|
| 280 | set_title = 'Static stability (K.m!U-1!N)' |
|---|
| 281 | set_titlex = 'Static stability (K.m!U-1!N)' |
|---|
| 282 | set_titley = 'Altitude above surface (km)' |
|---|
| 283 | set_subtitle = 'Mean over the simulation domain' |
|---|
| 284 | set_xrange = [-1.,5.] |
|---|
| 285 | set_yrange = [0.,10.] |
|---|
| 286 | set_tickx = 1. |
|---|
| 287 | set_ticky = 1. |
|---|
| 288 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 289 | localtimes = [9,11,13,15,17] |
|---|
| 290 | localtimes = [17] |
|---|
| 291 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 292 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 293 | !P.Charsize = 1.2 |
|---|
| 294 | altlin=0 & loadct, 0 |
|---|
| 295 | user_lt=localtimes[0] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 296 | 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 |
|---|
| 297 | oplot, zefield(*,nntt), zey(*,nntt), psym=5 |
|---|
| 298 | for ll = 1, n_elements(localtimes)-1 do begin |
|---|
| 299 | CASE altlin OF |
|---|
| 300 | 0: altlin=1 |
|---|
| 301 | 1: altlin=0 |
|---|
| 302 | ENDCASE |
|---|
| 303 | user_lt=localtimes[ll] & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) & nntt = yeah(0) |
|---|
| 304 | if (nntt ne -1) then oplot, zefield(*,nntt), zey(*,nntt), linestyle=altlin |
|---|
| 305 | ;if (nntt ne -1) then oplot, zefield(*,nntt), zey(*,nntt), psym=5 |
|---|
| 306 | endfor |
|---|
| 307 | oplot, 0.*zefield(*,nntt) + 1.5, zey(*,nntt), linestyle=2 |
|---|
| 308 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 309 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 310 | |
|---|
| 311 | |
|---|
| 312 | |
|---|
| 313 | stop |
|---|
| 314 | |
|---|
| 315 | |
|---|
| 316 | |
|---|
| 317 | |
|---|
| 318 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 319 | goto, no_staticstab |
|---|
| 320 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 321 | if (saveps eq 'true') then begin |
|---|
| 322 | device, /close |
|---|
| 323 | set_plot, 'ps' & device, filename='plot/staticstab.ps' |
|---|
| 324 | endif |
|---|
| 325 | |
|---|
| 326 | user_lt=17. & yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) |
|---|
| 327 | |
|---|
| 328 | ;;; recompute height @ given local time |
|---|
| 329 | caca=heightp+hgtu/1000. |
|---|
| 330 | height=reform(ph(*,0,*,*)) |
|---|
| 331 | height=total(height,1)/n_elements(height(*,0,0)) |
|---|
| 332 | height=reform(height) |
|---|
| 333 | height=reform(height(*,yeah(0))) |
|---|
| 334 | height=height/1000./3.72 |
|---|
| 335 | heightp=height(0:n_elements(height(*))-2) |
|---|
| 336 | print, 'new minus old', heightp - caca |
|---|
| 337 | ;;; recompute height @ given local time |
|---|
| 338 | |
|---|
| 339 | press=reform(p(*,0,*,*)) |
|---|
| 340 | press=total(press,1)/n_elements(press(*,0,0)) |
|---|
| 341 | press=reform(press) |
|---|
| 342 | press=reform(press(*,yeah(0))) |
|---|
| 343 | press=press(0:n_elements(press(*))-2) |
|---|
| 344 | |
|---|
| 345 | staticstab = DERIV( heightp , reform(what_I_plot5(*,yeah(0))) ) + 1000.*3.72/844.6 ;; 4.9 dans Hinson |
|---|
| 346 | staticstab = staticstab(0:n_elements(staticstab)-5) |
|---|
| 347 | heightp = heightp(0:n_elements(heightp)-5) |
|---|
| 348 | |
|---|
| 349 | plot, $ |
|---|
| 350 | staticstab, $ |
|---|
| 351 | heightp,$ |
|---|
| 352 | xtitle='Static stability (K/km)',$ |
|---|
| 353 | xrange=[-1.,5.], $ |
|---|
| 354 | xtickinterval=0.5, $ |
|---|
| 355 | ytitle='Altitude above surface (km)', $ |
|---|
| 356 | yrange=[min(heightp),max(heightp)], $ |
|---|
| 357 | ytickinterval=1., $ |
|---|
| 358 | title="LMD LES Static stability @ LT 17h (K/km)", $ |
|---|
| 359 | subtitle='zonal average at lat. '+latwrite |
|---|
| 360 | oplot, $ |
|---|
| 361 | staticstab, $ |
|---|
| 362 | heightp,$ |
|---|
| 363 | psym=5 |
|---|
| 364 | oplot, $ |
|---|
| 365 | findgen(n_elements(heightp))*0. + 1.,$ |
|---|
| 366 | heightp,$ |
|---|
| 367 | linestyle=2 |
|---|
| 368 | oplot, $ |
|---|
| 369 | findgen(n_elements(heightp))*0. + 2.,$ |
|---|
| 370 | heightp,$ |
|---|
| 371 | linestyle=2 |
|---|
| 372 | |
|---|
| 373 | ;;;;;;;;;; |
|---|
| 374 | w = where((staticstab gt 1.5) and ((heightp- hgtu/1000.) gt 1.)) |
|---|
| 375 | t_top_plus = what_I_plot5(w(0),yeah(0)) |
|---|
| 376 | z_top_plus = heightp(w(0)) |
|---|
| 377 | p_top_plus = press(w(0)) |
|---|
| 378 | t_top_moins = what_I_plot5(w(0)-1,yeah(0)) |
|---|
| 379 | z_top_moins = heightp(w(0)-1) |
|---|
| 380 | p_top_moins = press(w(0)-1) |
|---|
| 381 | 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. |
|---|
| 382 | xyouts, 3., 1.5 + (max(heightp) + min(heightp)) / 3., 'Ls = '+string(lsu,'(F5.1)')+'!Uo!N', CHARSIZE=1 |
|---|
| 383 | xyouts, 3., 1. + (max(heightp) + min(heightp)) / 3., 'Lat = '+string(latu,'(F5.1)')+'!Uo!N', CHARSIZE=1 |
|---|
| 384 | xyouts, 3., 0.5 + (max(heightp) + min(heightp)) / 3., 'LonE = '+string(lonu,'(F6.1)')+'!Uo!N', CHARSIZE=1 |
|---|
| 385 | xyouts, 3., (max(heightp) + min(heightp)) / 3., 'T!Dt!N = '+string(t_top_plus,'(I0)')+'/'+string(t_top_moins,'(I0)')+' K ', CHARSIZE=1 |
|---|
| 386 | 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 |
|---|
| 387 | 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 |
|---|
| 388 | xyouts, 3., -1.5 + (max(heightp) + min(heightp)) / 3., 'z!Ds!N = '+string(hgtu/1000.,'(F4.1)')+' km', CHARSIZE=1 |
|---|
| 389 | xyouts, 3., -2. + (max(heightp) + min(heightp)) / 3., 'D = '+string(pbl_depth,'(F4.1)')+' km', CHARSIZE=1 |
|---|
| 390 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 391 | no_staticstab: |
|---|
| 392 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 393 | |
|---|
| 394 | |
|---|
| 395 | |
|---|
| 396 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 397 | if (saveps eq 'true') then begin |
|---|
| 398 | device, /close |
|---|
| 399 | endif |
|---|
| 400 | stop |
|---|
| 401 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 402 | |
|---|
| 403 | |
|---|
| 404 | user_lt=13. |
|---|
| 405 | yeah=where(abs(localtime-user_lt) eq (min(abs(localtime-user_lt)))) |
|---|
| 406 | |
|---|
| 407 | user_h=0.1 |
|---|
| 408 | walt=where(abs(heightp-user_h) eq (min(abs(heightp-user_h)))) |
|---|
| 409 | |
|---|
| 410 | ;mapfield=reform(t[*,*,walt(0),w(0)]) |
|---|
| 411 | ;help, mapfield |
|---|
| 412 | ;contour, mapfield |
|---|
| 413 | |
|---|
| 414 | section=reform(w[*,0,*,yeah(0)]) |
|---|
| 415 | |
|---|
| 416 | section=reform(w[*,0,*,160]) |
|---|
| 417 | |
|---|
| 418 | lev=[-12.,-8.,-4.,4.,8.,12.] |
|---|
| 419 | lev=[-5.,-4.,-3.,-2.,-1.,1.,2.,3.,4.,5.] |
|---|
| 420 | contour, $ |
|---|
| 421 | section, $ |
|---|
| 422 | (lon(*,0)-lon(0,0))*59., $ |
|---|
| 423 | heightp, $ |
|---|
| 424 | levels=lev , $ |
|---|
| 425 | C_LINESTYLE = (lev LT 0.0) |
|---|
| 426 | |
|---|
| 427 | device, /close |
|---|
| 428 | |
|---|
| 429 | end |
|---|
| 430 | |
|---|
