| 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 | velmax = fltarr(nz,nttot) & velmin = fltarr(nz,nttot) |
|---|
| 15 | ustarm = fltarr(nttot) & ustarmax = fltarr(nttot) |
|---|
| 16 | for loop = 0, nf-1 do begin |
|---|
| 17 | timetime = SYSTIME(1) |
|---|
| 18 | if (loop ne nf-1) then nloop2=nt else nloop2=ntlast |
|---|
| 19 | for loop2 = 0, nloop2-1 do begin |
|---|
| 20 | uprime = getget(filesWRF(loop), 'U', count=[0,0,0,1], offset=[0,0,0,loop2]) |
|---|
| 21 | vprime = getget(filesWRF(loop), 'V', count=[0,0,0,1], offset=[0,0,0,loop2]) |
|---|
| 22 | ; |
|---|
| 23 | ; autre methode... ustar² = sqrt(<u'w'>s²+<v'w'>s²) |
|---|
| 24 | ; |
|---|
| 25 | !QUIET=1 |
|---|
| 26 | ustar = getget(filesWRF(loop), 'USTM', count=[0,0,0,1], offset=[0,0,0,loop2]) |
|---|
| 27 | ustarmax(yeye) = max(reform(ustar(*,*,0,0))) |
|---|
| 28 | ustarm(yeye) = reform(TOTAL(TOTAL(TEMPORARY(ustar),1),1) / float(nx) / float(ny)) |
|---|
| 29 | print, ustarm(yeye), ustarmax(yeye) |
|---|
| 30 | !QUIET=0 |
|---|
| 31 | veltot = sqrt(uprime^2+vprime^2) ;; attention si background wind |
|---|
| 32 | for toto = 0, nz-1 do begin |
|---|
| 33 | velmax(toto,yeye) = max(reform(veltot(*,*,toto,0))) |
|---|
| 34 | ;print, velmax(toto,yeye) |
|---|
| 35 | endfor |
|---|
| 36 | yeye = TEMPORARY(yeye) + 1 |
|---|
| 37 | endfor |
|---|
| 38 | endfor |
|---|
| 39 | save, ustarm, ustarmax, velmax, localtime, filename='addturb.dat' ;; pour l'instant ecrase systematiquement |
|---|
| 40 | stop |
|---|
| 41 | ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
|---|
| 42 | |
|---|
