| 1 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 2 | OPENR, 22, 'input_coord' & READF, 22, lonu & READF, 22, latu & READF, 22, lsu & READF, 22, lctu & CLOSE, 22 |
|---|
| 3 | OPENR, 23, 'input_more' & READF, 23, hgtu, tsurfu & CLOSE, 23 |
|---|
| 4 | domain='d01' & filesWRF = FindFile('wrfout_'+domain+'_????-??-??_??:??:??') & nf=n_elements(filesWRF) |
|---|
| 5 | id=ncdf_open(filesWRF(0)) |
|---|
| 6 | NCDF_DIMINQ, id, NCDF_DIMID(id, 'west_east' ), toto, nx & NCDF_DIMINQ, id, NCDF_DIMID(id, 'south_north' ), toto, ny |
|---|
| 7 | NCDF_DIMINQ, id, NCDF_DIMID(id, 'bottom_top' ), toto, nz & NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), toto, nt |
|---|
| 8 | NCDF_CLOSE, id |
|---|
| 9 | id=ncdf_open(filesWRF(nf-1)) ;; for interrupted runs |
|---|
| 10 | NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time' ), toto, ntlast |
|---|
| 11 | NCDF_CLOSE, id |
|---|
| 12 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 13 | yeye = 0 & nttot = (nf-1)*nt + ntlast & localtime = lctu + history_interval_s*findgen(nttot)/3700. |
|---|
| 14 | std_tprime = fltarr(nz,nttot) |
|---|
| 15 | field = fltarr(nz,nttot) |
|---|
| 16 | anomalt = 1. |
|---|
| 17 | for loop = 0, nf-1 do begin |
|---|
| 18 | timetime = SYSTIME(1) |
|---|
| 19 | print, loop |
|---|
| 20 | if (loop ne nf-1) then nloop2=nt else nloop2=ntlast |
|---|
| 21 | for loop2 = 0, nloop2-1 do begin |
|---|
| 22 | !QUIET=1 |
|---|
| 23 | |
|---|
| 24 | ;;;; TEMPPOT |
|---|
| 25 | ;tprime = getget(filesWRF(loop), 'T', anomaly=anomalt, count=[0,0,0,1], offset=[0,0,0,loop2]) ;; t' = t - <t> |
|---|
| 26 | ;;;; TEMP (K) |
|---|
| 27 | tempk = reform ( ( t0 + getget(filesWRF(loop), 'T', count=[0,0,0,1], offset=[0,0,0,loop2]) ) * ( getget(filesWRF(loop), 'PTOT', count=[0,0,0,1], offset=[0,0,0,loop2]) / p0 )^r_cp ) |
|---|
| 28 | meanmean = reform(TOTAL(TOTAL(tempk,1),1) / float(nx) / float(ny)) |
|---|
| 29 | for i=0,n_elements(tempk(*,0,0,0))-1 do for j=0,n_elements(tempk(0,*,0,0))-1 do tempk(i,j,*) = tempk(i,j,*) - meanmean |
|---|
| 30 | tprime = TEMPORARY(tempk)^2 |
|---|
| 31 | |
|---|
| 32 | std_tprime(*, yeye) = sqrt ( reform(TOTAL(TOTAL(TEMPORARY(tprime),1),1) / float(nx) / float(ny)) ) |
|---|
| 33 | ;print, std_tprime(*, yeye) |
|---|
| 34 | ;print, yeye |
|---|
| 35 | |
|---|
| 36 | !QUIET=0 |
|---|
| 37 | yeye = TEMPORARY(yeye) + 1 |
|---|
| 38 | endfor |
|---|
| 39 | endfor |
|---|
| 40 | ;save, std_tprime, localtime, filename='addturb.dat' ;; pour l'instant ecrase systematiquement |
|---|
| 41 | ;stop |
|---|
| 42 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 43 | |
|---|
| 44 | |
|---|
| 45 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 46 | zefield = transpose(std_tprime) |
|---|
| 47 | zex = localtime |
|---|
| 48 | zey = h |
|---|
| 49 | set_name = 'std_tprime.ps' |
|---|
| 50 | set_title = "Temperature standard deviation (K)" ;"Maximum updraft speed (m.s!U-1!N)" |
|---|
| 51 | set_titlex = 'Local Time (h)' |
|---|
| 52 | set_titley = 'Altitude above surface (km)' |
|---|
| 53 | set_subtitle = '' |
|---|
| 54 | set_xrange = [8.,18.] |
|---|
| 55 | set_yrange = [0.,7.] |
|---|
| 56 | set_tickx = 1. |
|---|
| 57 | set_ticky = 1. |
|---|
| 58 | minval = 0. |
|---|
| 59 | maxval = 2. ;15.;14.;13. |
|---|
| 60 | nlev = 10 ;(maxval-minval)*10. |
|---|
| 61 | pal = 22 |
|---|
| 62 | rrr = 'no' |
|---|
| 63 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 64 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 65 | !P.Charsize = 1.2 |
|---|
| 66 | ;; 0. levels |
|---|
| 67 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 68 | ;;; 1. background |
|---|
| 69 | 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 |
|---|
| 70 | ;; 2. color field |
|---|
| 71 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 72 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 73 | ;; 3. contour field |
|---|
| 74 | 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) |
|---|
| 75 | ;; 4. choose output |
|---|
| 76 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 77 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 78 | |
|---|
| 79 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 80 | zefield = transpose(std_tprime) |
|---|
| 81 | zex = localtime |
|---|
| 82 | zey = h |
|---|
| 83 | set_name = 'std_tprime_ns.ps' |
|---|
| 84 | set_title = "Temperature standard deviation (K)" ;"Maximum updraft speed (m.s!U-1!N)" |
|---|
| 85 | set_titlex = 'Local Time (h)' |
|---|
| 86 | set_titley = 'Altitude above surface (km)' |
|---|
| 87 | set_subtitle = '' |
|---|
| 88 | set_xrange = [8.,18.] |
|---|
| 89 | set_yrange = [0.,1.] |
|---|
| 90 | set_tickx = 1. |
|---|
| 91 | set_ticky = 1. |
|---|
| 92 | minval = 0. |
|---|
| 93 | maxval = 2. ;15.;14.;13. |
|---|
| 94 | nlev = 10 ;(maxval-minval)*10. |
|---|
| 95 | pal = 22 |
|---|
| 96 | rrr = 'no' |
|---|
| 97 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 98 | if (saveps eq 'true') then PS_Start, FILENAME=set_name |
|---|
| 99 | !P.Charsize = 1.2 |
|---|
| 100 | ;; 0. levels |
|---|
| 101 | lev = minval + (maxval-minval)*findgen(nlev+1)/float(nlev) & if (minval ne 0.) then lev = lev[where(lev ne 0.)] |
|---|
| 102 | ;;; 1. background |
|---|
| 103 | 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 |
|---|
| 104 | ;; 2. color field |
|---|
| 105 | loadct, pal & if (rrr eq 'yes') then TVLCT, r, g, b, /Get & if (rrr eq 'yes') then TVLCT, Reverse(r), Reverse(g), Reverse(b) |
|---|
| 106 | contour, zefield, zex, zey, levels=lev, c_labels=findgen(n_elements(lev))*0.+1., /overplot, /cell_fill |
|---|
| 107 | ;; 3. contour field |
|---|
| 108 | 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) |
|---|
| 109 | ;; 4. choose output |
|---|
| 110 | if (saveps eq 'true') then PS_End, /PNG |
|---|
| 111 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 112 | |
|---|
| 113 | |
|---|
| 114 | |
|---|
| 115 | |
|---|
| 116 | stop |
|---|
