source: trunk/LMDZ.TITAN/Tools/fft.F90 @ 819

Last change on this file since 819 was 816, checked in by slebonnois, 12 years ago

SL: tools for postprocessing (Veznus and Titan); see DOC/documentation/vt-tools.pdf
File size: 28.2 KB
Line 
1program fft
2
3! SL 01/2010:
4! This program reads 4D (lon-lat-alt-time) fields recast in log P coordinates
5!
6! it computes fft of temperature, zonal and merid winds from high-frequency outputs:
7!
8! fftaT -- 4D -- FFT in amplitude of temperature field (K)
9! fftau -- 4D -- FFT in amplitude of zonal wind (m s-1)
10! fftav -- 4D -- FFT in amplitude of meridional wind (m s-1)
11! ulf   -- 4D -- low  freq part of zonal wind perturbation uprim (m s-1)
12! ubf   -- 4D -- band freq part of zonal wind perturbation uprim (m s-1)
13! uhf   -- 4D -- high freq part of zonal wind perturbation uprim (m s-1)
14! vlf   -- 4D -- low  freq part of meridional wind perturbation vprim (m s-1)
15! vbf   -- 4D -- band freq part of meridional wind perturbation vprim (m s-1)
16! vhf   -- 4D -- high freq part of meridional wind perturbation vprim (m s-1)
17! wlf   -- 4D -- low  freq part of vertical wind perturbation wprim (Pa s-1)
18! wbf   -- 4D -- band freq part of vertical wind perturbation wprim (Pa s-1)
19! whf   -- 4D -- high freq part of vertical wind perturbation wprim (Pa s-1)
20! Tlf   -- 4D -- low  freq part of temperature perturbation Tprim (K)
21! Tbf   -- 4D -- band freq part of temperature perturbation Tprim (K)
22! Thf   -- 4D -- high freq part of temperature perturbation Tprim (K)
23!
24! Minimal requirements and dependencies:
25! The dataset must include the following data:
26! - pressure vertical coordinate
27! - atmospheric temperature
28! - zonal, meridional and vertical winds
29!
30! We use the FFTW library:   http://www.fftw.org
31! These routines are in C, but also include Fortran interfaces.
32!
33! Convention: qbar  <=> zonal average    / qstar = q - qbar
34!             qmean <=> temporal average / qprim = q - qmean
35!
36!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
37!  FILTRES
38!!!!!!!!!!!!!!!!!!!!!!!!!!!
39!  low  frequencies: qlf= lowfreq(qprim)
40!  band frequencies: qbf=bandfreq(qprim)
41!  high frequencies: qhf=highfreq(qprim)
42!
43!  Les frequences seuils sont ajustables dans filter.h
44!
45!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
46
47implicit none
48
49include "netcdf.inc" ! NetCDF definitions
50
51character (len=128) :: infile ! input file name (name_P.nc)
52character (len=128) :: outfile1,outfile2,outfile3,outfile4 ! output file names
53
54character (len=64) :: text ! to store some text
55integer infid ! NetCDF input file ID
56integer outfid1,outfid2,outfid3,outfid4 ! NetCDF output files ID
57integer lon_dimid1,lat_dimid1,alt_dimid1,time_dimid1 ! NetCDF dimension IDs
58integer lon_dimid2,lat_dimid2,alt_dimid2,time_dimid2 ! NetCDF dimension IDs
59integer lon_dimid3,lat_dimid3,alt_dimid3,time_dimid3 ! NetCDF dimension IDs
60integer lon_dimid4,lat_dimid4,alt_dimid4,time_dimid4 ! NetCDF dimension IDs
61integer lon_varid,lat_varid,alt_varid,time_varid
62integer              :: datashape1d ! shape of 1D datasets
63integer,dimension(4) :: datashape4d ! shape of 4D datasets
64
65real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
66real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
67real :: pi
68real,dimension(:),allocatable :: lon ! longitude
69integer lonlength ! # of grid points along longitude
70real,dimension(:),allocatable :: lat ! latitude
71real,dimension(:),allocatable :: latrad ! latitude in radian
72integer latlength ! # of grid points along latitude
73real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
74integer altlength ! # of grid point along altitude (of input datasets)
75real,dimension(:),allocatable :: time ! time
76real,dimension(:),allocatable :: freq ! frequencies of the FFT (only timelength/2+1 values)
77integer timelength ! # of points along time
78real,dimension(:,:,:,:),allocatable :: temp ! atmospheric temperature
79real,dimension(:,:,:,:),allocatable :: vitu ! zonal wind (in m/s)
80real,dimension(:,:,:,:),allocatable :: vitv ! meridional wind (in m/s)
81real,dimension(:,:,:,:),allocatable :: vitw ! vertical wind (in Pa/s)
82
83!!! output variables
84real,dimension(:,:,:,:),allocatable :: fftaT ! FFT in amplitude of temperature (K)
85real,dimension(:,:,:,:),allocatable :: fftau ! FFT in amplitude of zonal wind (m s-1)
86real,dimension(:,:,:,:),allocatable :: fftav ! FFT in amplitude of meridional wind (m s-1)
87real,dimension(:,:,:,:),allocatable :: fftaw ! FFT in amplitude of vertical wind (Pa s-1)
88real,dimension(:,:,:,:),allocatable :: ulf ! low  freq part of zonal wind perturbation uprim (m s-1)
89real,dimension(:,:,:,:),allocatable :: ubf ! band freq part of zonal wind perturbation uprim (m s-1)
90real,dimension(:,:,:,:),allocatable :: uhf ! high freq part of zonal wind perturbation uprim (m s-1)
91real,dimension(:,:,:,:),allocatable :: vlf ! low  freq part of meridional wind perturbation vprim (m s-1)
92real,dimension(:,:,:,:),allocatable :: vbf ! band freq part of meridional wind perturbation vprim (m s-1)
93real,dimension(:,:,:,:),allocatable :: vhf ! high freq part of meridional wind perturbation vprim (m s-1)
94real,dimension(:,:,:,:),allocatable :: wlf ! low  freq part of vertical wind perturbation vprim (Pa s-1)
95real,dimension(:,:,:,:),allocatable :: wbf ! band freq part of vertical wind perturbation vprim (Pa s-1)
96real,dimension(:,:,:,:),allocatable :: whf ! high freq part of vertical wind perturbation vprim (Pa s-1)
97real,dimension(:,:,:,:),allocatable :: Tlf ! low  freq part of temperature perturbation Tprim (K)
98real,dimension(:,:,:,:),allocatable :: Tbf ! band freq part of temperature perturbation Tprim (K)
99real,dimension(:,:,:,:),allocatable :: Thf ! high freq part of temperature perturbation Tprim (K)
100
101! local variables
102real,dimension(:,:,:,:),allocatable :: uprim
103real,dimension(:,:,:,:),allocatable :: vprim
104real,dimension(:,:,:,:),allocatable :: wprim
105real,dimension(:,:,:,:),allocatable :: Tprim
106
107! lon,lat,alt
108real,dimension(:,:,:),allocatable :: umean
109real,dimension(:,:,:),allocatable :: vmean
110real,dimension(:,:,:),allocatable :: wmean
111real,dimension(:,:,:),allocatable :: Tmean
112
113! for FFTW routines
114real,dimension(:),allocatable :: wndow
115double precision,dimension(:),allocatable :: var,fltvar
116double complex,dimension(:),allocatable :: fftvar,fltfft
117double complex,dimension(:),allocatable :: filtrelf,filtrebf,filtrehf
118integer   :: M_fft
119integer*8 :: planf,planb
120
121
122integer ierr,ierr1,ierr2 ! NetCDF routines return codes
123integer i,j,ilon,ilat,ilev,itim ! for loops
124logical flagfft
125logical :: lmdflag
126
127include "planet.h"
128include "filter.h"
129
130#include <fftw3.f>
131
132!===============================================================================
133! 1. Input parameters
134!===============================================================================
135
136pi = 2.*asin(1.)
137
138write(*,*) ""
139write(*,*) "You are working on the atmosphere of ",planet
140
141!===============================================================================
142! 1.1 Input file
143!===============================================================================
144
145write(*,*) ""
146write(*,*) " Program valid for files with pressure axis (*_P.nc)"
147write(*,*) "Enter input file name:"
148
149read(*,'(a128)') infile
150write(*,*) ""
151
152! open input file
153
154ierr = NF_OPEN(infile,NF_NOWRITE,infid)
155if (ierr.ne.NF_NOERR) then
156   write(*,*) 'ERROR: Pb opening file ',trim(infile)
157   stop ""
158endif
159
160!===============================================================================
161! 1.2 Get grids in lon,lat,alt(pressure),time
162!===============================================================================
163
164call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
165                     alt_varid,altlength,time_varid,timelength,lmdflag )
166
167allocate(lon(lonlength))
168ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
169if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
170
171allocate(lat(latlength))
172ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
173if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
174
175allocate(latrad(latlength))
176latrad = lat*pi/180.
177
178allocate(plev(altlength))
179ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
180if (ierr.ne.NF_NOERR) stop "Error: Failed reading altitude (ie pressure levels)"
181
182allocate(time(timelength))
183ierr=NF_GET_VAR_REAL(infid,time_varid,time)
184if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
185
186!===============================================================================
187! 1.3 Get output file name
188!===============================================================================
189write(*,*) ""
190!write(*,*) "Enter output file name"
191!read(*,*) outfile
192outfile1=infile(1:len_trim(infile)-3)//"_UFFT.nc"
193outfile2=infile(1:len_trim(infile)-3)//"_VFFT.nc"
194outfile3=infile(1:len_trim(infile)-3)//"_WFFT.nc"
195outfile4=infile(1:len_trim(infile)-3)//"_TFFT.nc"
196write(*,*) "Output file names are: "
197if (ok_out(1)) write(*,*) trim(outfile1)
198if (ok_out(2)) write(*,*) trim(outfile2)
199if (ok_out(3)) write(*,*) trim(outfile3)
200if (ok_out(4)) write(*,*) trim(outfile4)
201
202
203!===============================================================================
204! 2.1 Store needed fields
205!===============================================================================
206
207!===============================================================================
208! 2.1.1 Atmospheric temperature
209!===============================================================================
210allocate(temp(lonlength,latlength,altlength,timelength))
211
212text="temp"
213call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
214if (ierr1.ne.NF_NOERR) then
215  write(*,*) "  looking for t instead... "
216  text="t"
217  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
218  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
219endif
220if (ierr2.ne.NF_NOERR) stop "Error: Failed reading temperature"
221
222!===============================================================================
223! 2.1.2 Winds
224!===============================================================================
225allocate(vitu(lonlength,latlength,altlength,timelength))
226allocate(vitv(lonlength,latlength,altlength,timelength))
227allocate(vitw(lonlength,latlength,altlength,timelength))
228
229! zonal wind vitu (in m/s)
230text="vitu"
231call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
232if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
233if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
234
235! meridional wind vitv (in m/s)
236text="vitv"
237call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
238if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
239if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
240
241! vertical wind vitw (in Pa/s)
242text="vitw"
243call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
244if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
245if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
246
247!===============================================================================
248! 2.2 Computations
249!===============================================================================
250
251print*,"debut calcul"
252
253!===============================================================================
254! 2.2.1 FFT and filtering
255!===============================================================================
256! allocations
257!-------------
258allocate(fftau(lonlength,latlength,altlength,timelength))
259allocate(fftav(lonlength,latlength,altlength,timelength))
260allocate(fftaw(lonlength,latlength,altlength,timelength))
261allocate(fftaT(lonlength,latlength,altlength,timelength))
262allocate(uprim(lonlength,latlength,altlength,timelength))
263allocate(vprim(lonlength,latlength,altlength,timelength))
264allocate(wprim(lonlength,latlength,altlength,timelength))
265allocate(Tprim(lonlength,latlength,altlength,timelength))
266allocate(ulf(lonlength,latlength,altlength,timelength))
267allocate(vlf(lonlength,latlength,altlength,timelength))
268allocate(wlf(lonlength,latlength,altlength,timelength))
269allocate(Tlf(lonlength,latlength,altlength,timelength))
270allocate(ubf(lonlength,latlength,altlength,timelength))
271allocate(vbf(lonlength,latlength,altlength,timelength))
272allocate(wbf(lonlength,latlength,altlength,timelength))
273allocate(Tbf(lonlength,latlength,altlength,timelength))
274allocate(uhf(lonlength,latlength,altlength,timelength))
275allocate(vhf(lonlength,latlength,altlength,timelength))
276allocate(whf(lonlength,latlength,altlength,timelength))
277allocate(Thf(lonlength,latlength,altlength,timelength))
278
279! lon,lat,alt
280allocate(umean(lonlength,latlength,altlength))
281allocate(vmean(lonlength,latlength,altlength))
282allocate(wmean(lonlength,latlength,altlength))
283allocate(Tmean(lonlength,latlength,altlength))
284
285! time / frequencies
286allocate(freq(timelength))
287allocate(wndow(timelength))
288allocate(var(timelength))
289allocate(fltvar(timelength))
290M_fft = timelength/2
291allocate(fftvar(M_fft+1))
292allocate(fltfft(M_fft+1))
293allocate(filtrelf(M_fft+1))
294allocate(filtrebf(M_fft+1))
295allocate(filtrehf(M_fft+1))
296
297! intermediates
298!-----------------
299
300call moytim(lonlength,latlength,altlength,timelength,miss_val,vitu,umean)
301call moytim(lonlength,latlength,altlength,timelength,miss_val,vitv,vmean)
302call moytim(lonlength,latlength,altlength,timelength,miss_val,vitw,wmean)
303call moytim(lonlength,latlength,altlength,timelength,miss_val,temp,Tmean)
304
305do ilon=1,lonlength
306 do ilat=1,latlength
307  do ilev=1,altlength
308   do itim=1,timelength
309    if ((vitu(ilon,ilat,ilev,itim).ne.miss_val).and. &
310        (umean(ilon,ilat,ilev)    .ne.miss_val)) then
311  uprim(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim)-umean(ilon,ilat,ilev)
312    else
313  uprim(ilon,ilat,ilev,itim) = miss_val
314    endif
315    if ((vitv(ilon,ilat,ilev,itim).ne.miss_val).and. &
316        (vmean(ilon,ilat,ilev)    .ne.miss_val)) then
317  vprim(ilon,ilat,ilev,itim) = vitv(ilon,ilat,ilev,itim)-vmean(ilon,ilat,ilev)
318    else
319  vprim(ilon,ilat,ilev,itim) = miss_val
320    endif
321    if ((vitw(ilon,ilat,ilev,itim).ne.miss_val).and. &
322        (wmean(ilon,ilat,ilev)    .ne.miss_val)) then
323  wprim(ilon,ilat,ilev,itim) = vitw(ilon,ilat,ilev,itim)-wmean(ilon,ilat,ilev)
324    else
325  wprim(ilon,ilat,ilev,itim) = miss_val
326    endif
327    if ((temp(ilon,ilat,ilev,itim).ne.miss_val).and. &
328        (Tmean(ilon,ilat,ilev)    .ne.miss_val)) then
329  Tprim(ilon,ilat,ilev,itim) = temp(ilon,ilat,ilev,itim)-Tmean(ilon,ilat,ilev)
330    else
331  Tprim(ilon,ilat,ilev,itim) = miss_val
332    endif
333   enddo
334  enddo
335 enddo
336enddo ! lonlength
337
338! fft intermediates
339!-------------
340
341! Define the frequencies
342do itim=1,M_fft+1
343  freq(itim) = (itim-1)/(timelength*(time(2)-time(1)))
344enddo
345do itim=M_fft+2,timelength
346  freq(itim) = 0.
347enddo
348
349! Define the window (triangle)
350do itim=1,timelength
351! N window:
352!  wndow(itim)= 1.
353! triangulaire de moyenne = 1
354  wndow(itim)= 2.*(1. - abs(real(itim-0.5-M_fft)/real(M_fft)))
355enddo
356
357! Define the filters
358
359! IN filter.h :
360! Low  cutting frequency, in Hz    : fcoup1
361! High cutting frequency, in Hz    : fcoup2
362! Half-width of the filters, in Hz : width
363
364!print*,"Low  cutting frequency, in Hz ?"
365!read(*,*) fcoup1
366!print*,"High cutting frequency, in Hz ?"
367!read(*,*) fcoup2
368!print*,"Half-width of the filters, in Hz ?"
369!read(*,*) width
370!width = 3./time(timelength)
371
372do itim=1,M_fft+1
373  if (freq(itim).lt.(fcoup1-width)) then
374     filtrelf(itim) = 1.
375  elseif (freq(itim).gt.(fcoup1+width)) then
376     filtrelf(itim) = 0.
377  else
378     filtrelf(itim) = (1.+sin(pi*(fcoup1-freq(itim))/(2.*width)))/2.
379  endif
380  if (freq(itim).lt.(fcoup2-width)) then
381     filtrehf(itim) = 0.
382  elseif (freq(itim).gt.(fcoup2+width)) then
383     filtrehf(itim) = 1.
384  else
385     filtrehf(itim) = (1.-sin(pi*(fcoup2-freq(itim))/(2.*width)))/2.
386  endif
387  filtrebf(itim) = (1.-filtrelf(itim))*(1.-filtrehf(itim))
388enddo
389
390
391! fft and filtering
392!-------------
393
394!---FFTW routines
395call dfftw_plan_dft_r2c_1d(planf,timelength,var,fftvar,FFTW_MEASURE)
396call dfftw_plan_dft_c2r_1d(planb,timelength,fltfft,fltvar,FFTW_MEASURE)
397!---
398
399do ilon=1,lonlength
400 do ilat=1,latlength
401  do ilev=1,altlength
402
403! For zonal wind field
404if (ok_out(1)) then
405
406   flagfft=.true.
407   do itim=1,timelength
408     if (uprim(ilon,ilat,ilev,itim).eq.miss_val) flagfft=.false.
409   enddo
410
411   if (flagfft) then
412
413! 1/ windowing to improve spectral analysis
414      var(:)=uprim(ilon,ilat,ilev,:)*wndow(:)
415! 2/ FFT computation
416!---FFTW routines
417      call dfftw_execute_dft_r2c(planf,var,fftvar)
418!---
419! 3/ Amplitude of the FFT, for spectral analysis
420      fftau(ilon,ilat,ilev,1)=abs(fftvar(1))/M_fft
421      do itim=2,M_fft
422       fftau(ilon,ilat,ilev,itim) = abs(fftvar(itim))/M_fft
423      enddo
424      fftau(ilon,ilat,ilev,M_fft+1)=abs(fftvar(M_fft+1))/M_fft
425      do itim=M_fft+2,timelength
426       fftau(ilon,ilat,ilev,itim) = 0.
427      enddo
428
429! 4/ filtering the FFT in three regions
430! filtering + normalisation (low freq)
431      fltfft(:) = fftvar(:)*filtrelf(:)/timelength
432! 5/ backward FFT for each region
433!---FFTW routines
434      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
435!---
436! 6/ reverse the windowing
437      ulf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
438
439! filtering + normalisation (band freq)
440      fltfft(:) = fftvar(:)*filtrebf(:)/timelength
441!---FFTW routines
442      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
443!---
444      ubf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
445
446! filtering + normalisation (high freq)
447      fltfft(:) = fftvar(:)*filtrehf(:)/timelength
448!---FFTW routines
449      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
450!---
451      uhf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
452
453   else
454     fftau(ilon,ilat,ilev,itim) = miss_val
455       ulf(ilon,ilat,ilev,itim) = miss_val
456       ubf(ilon,ilat,ilev,itim) = miss_val
457       uhf(ilon,ilat,ilev,itim) = miss_val
458   endif ! flagfft
459
460endif !ok_out(1)
461
462! For meridional wind wind field
463if (ok_out(2)) then
464
465   flagfft=.true.
466   do itim=1,timelength
467     if (vprim(ilon,ilat,ilev,itim).eq.miss_val) flagfft=.false.
468   enddo
469
470   if (flagfft) then
471
472! 1/ windowing to improve spectral analysis
473      var(:)=vprim(ilon,ilat,ilev,:)*wndow(:)
474! 2/ FFT computation
475!---FFTW routines
476      call dfftw_execute_dft_r2c(planf,var,fftvar)
477!---
478! 3/ Amplitude of the FFT, for spectral analysis
479      fftav(ilon,ilat,ilev,1)=abs(fftvar(1))/M_fft
480      do itim=2,M_fft
481       fftav(ilon,ilat,ilev,itim) = abs(fftvar(itim))/M_fft
482      enddo
483      fftav(ilon,ilat,ilev,M_fft+1)=abs(fftvar(M_fft+1))/M_fft
484      do itim=M_fft+2,timelength
485       fftav(ilon,ilat,ilev,itim) = 0.
486      enddo
487
488! 4/ filtering the FFT in three regions
489! filtering + normalisation (low freq)
490      fltfft(:) = fftvar(:)*filtrelf(:)/timelength
491! 5/ backward FFT for each region
492!---FFTW routines
493      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
494!---
495! 6/ reverse the windowing
496      vlf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
497
498! filtering + normalisation (band freq)
499      fltfft(:) = fftvar(:)*filtrebf(:)/timelength
500!---FFTW routines
501      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
502!---
503      vbf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
504
505! filtering + normalisation (high freq)
506      fltfft(:) = fftvar(:)*filtrehf(:)/timelength
507!---FFTW routines
508      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
509!---
510      vhf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
511
512   else
513     fftav(ilon,ilat,ilev,itim) = miss_val
514       vlf(ilon,ilat,ilev,itim) = miss_val
515       vbf(ilon,ilat,ilev,itim) = miss_val
516       vhf(ilon,ilat,ilev,itim) = miss_val
517   endif ! flagfft
518
519endif !ok_out(2)
520
521! For vertical wind wind field
522if (ok_out(3)) then
523
524   flagfft=.true.
525   do itim=1,timelength
526     if (wprim(ilon,ilat,ilev,itim).eq.miss_val) flagfft=.false.
527   enddo
528
529   if (flagfft) then
530
531! 1/ windowing to improve spectral analysis
532      var(:)=wprim(ilon,ilat,ilev,:)*wndow(:)
533! 2/ FFT computation
534!---FFTW routines
535      call dfftw_execute_dft_r2c(planf,var,fftvar)
536!---
537! 3/ Amplitude of the FFT, for spectral analysis
538      fftaw(ilon,ilat,ilev,1)=abs(fftvar(1))/M_fft
539      do itim=2,M_fft
540       fftaw(ilon,ilat,ilev,itim) = abs(fftvar(itim))/M_fft
541      enddo
542      fftaw(ilon,ilat,ilev,M_fft+1)=abs(fftvar(M_fft+1))/M_fft
543      do itim=M_fft+2,timelength
544       fftaw(ilon,ilat,ilev,itim) = 0.
545      enddo
546
547! 4/ filtering the FFT in three regions
548! filtering + normalisation (low freq)
549      fltfft(:) = fftvar(:)*filtrelf(:)/timelength
550! 5/ backward FFT for each region
551!---FFTW routines
552      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
553!---
554! 6/ reverse the windowing
555      wlf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
556
557! filtering + normalisation (band freq)
558      fltfft(:) = fftvar(:)*filtrebf(:)/timelength
559!---FFTW routines
560      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
561!---
562      wbf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
563
564! filtering + normalisation (high freq)
565      fltfft(:) = fftvar(:)*filtrehf(:)/timelength
566!---FFTW routines
567      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
568!---
569      whf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
570
571   else
572     fftaw(ilon,ilat,ilev,itim) = miss_val
573       wlf(ilon,ilat,ilev,itim) = miss_val
574       wbf(ilon,ilat,ilev,itim) = miss_val
575       whf(ilon,ilat,ilev,itim) = miss_val
576   endif ! flagfft
577
578endif !ok_out(3)
579
580! For temperature field
581if (ok_out(4)) then
582
583   flagfft=.true.
584   do itim=1,timelength
585     if (Tprim(ilon,ilat,ilev,itim).eq.miss_val) flagfft=.false.
586   enddo
587
588   if (flagfft) then
589
590! 1/ windowing to improve spectral analysis
591      var(:)=Tprim(ilon,ilat,ilev,:)*wndow(:)
592! 2/ FFT computation
593!---FFTW routines
594      call dfftw_execute_dft_r2c(planf,var,fftvar)
595!---
596! 3/ Amplitude of the FFT, for spectral analysis
597      fftaT(ilon,ilat,ilev,1)=abs(fftvar(1))/M_fft
598      do itim=2,M_fft
599       fftaT(ilon,ilat,ilev,itim) = abs(fftvar(itim))/M_fft
600      enddo
601      fftaT(ilon,ilat,ilev,M_fft+1)=abs(fftvar(M_fft+1))/M_fft
602      do itim=M_fft+2,timelength
603       fftaT(ilon,ilat,ilev,itim) = 0.
604      enddo
605
606! 4/ filtering the FFT in three regions
607! filtering + normalisation (low freq)
608      fltfft(:) = fftvar(:)*filtrelf(:)/timelength
609! 5/ backward FFT for each region
610!---FFTW routines
611      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
612!---
613! 6/ reverse the windowing
614      Tlf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
615
616! filtering + normalisation (band freq)
617      fltfft(:) = fftvar(:)*filtrebf(:)/timelength
618!---FFTW routines
619      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
620!---
621      Tbf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
622
623! filtering + normalisation (high freq)
624      fltfft(:) = fftvar(:)*filtrehf(:)/timelength
625!---FFTW routines
626      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
627!---
628      Thf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
629
630   else
631     fftaT(ilon,ilat,ilev,itim) = miss_val
632       Tlf(ilon,ilat,ilev,itim) = miss_val
633       Tbf(ilon,ilat,ilev,itim) = miss_val
634       Thf(ilon,ilat,ilev,itim) = miss_val
635   endif ! flagfft
636
637endif !ok_out(4)
638
639  enddo
640 enddo
641enddo ! lonlength
642
643!---FFTW routines
644call dfftw_destroy_plan(planf)
645call dfftw_destroy_plan(planb)
646!---
647
648print*,"End of computations"
649
650!===============================================================================
651! 3. Create output files
652!===============================================================================
653
654! Create output files
655if (ok_out(1)) then
656ierr=NF_CREATE(outfile1,NF_CLOBBER,outfid1)
657if (ierr.ne.NF_NOERR) then
658  write(*,*)"Error: could not create file ",outfile1
659  stop
660endif
661endif !ok_out(1)
662
663if (ok_out(2)) then
664ierr=NF_CREATE(outfile2,NF_CLOBBER,outfid2)
665if (ierr.ne.NF_NOERR) then
666  write(*,*)"Error: could not create file ",outfile2
667  stop
668endif
669endif !ok_out(2)
670
671if (ok_out(3)) then
672ierr=NF_CREATE(outfile3,NF_CLOBBER,outfid3)
673if (ierr.ne.NF_NOERR) then
674  write(*,*)"Error: could not create file ",outfile3
675  stop
676endif
677endif !ok_out(3)
678
679if (ok_out(4)) then
680ierr=NF_CREATE(outfile4,NF_CLOBBER,outfid4)
681if (ierr.ne.NF_NOERR) then
682  write(*,*)"Error: could not create file ",outfile4
683  stop
684endif
685endif !ok_out(4)
686
687!===============================================================================
688! 3.1. Define and write dimensions
689!===============================================================================
690
691if (ok_out(1)) &
692call write_dim(outfid1,lonlength,latlength,altlength,timelength, &
693    lon,lat,plev,time,lon_dimid1,lat_dimid1,alt_dimid1,time_dimid1)
694if (ok_out(2)) &
695call write_dim(outfid2,lonlength,latlength,altlength,timelength, &
696    lon,lat,plev,time,lon_dimid2,lat_dimid2,alt_dimid2,time_dimid2)
697if (ok_out(3)) &
698call write_dim(outfid3,lonlength,latlength,altlength,timelength, &
699    lon,lat,plev,time,lon_dimid3,lat_dimid3,alt_dimid3,time_dimid3)
700if (ok_out(4)) &
701call write_dim(outfid4,lonlength,latlength,altlength,timelength, &
702    lon,lat,plev,time,lon_dimid4,lat_dimid4,alt_dimid4,time_dimid4)
703
704!===============================================================================
705! 3.2. Define and write variables
706!===============================================================================
707
708if (ok_out(1)) then
709
710datashape4d(1)=lon_dimid1
711datashape4d(2)=lat_dimid1
712datashape4d(3)=alt_dimid1
713datashape4d(4)=time_dimid1
714datashape1d   =time_dimid1
715
716call write_var1d(outfid1,datashape1d,timelength,&
717                "freq      ", "FFT frequencies     ","s-1       ",miss_val,&
718                 freq )
719
720call write_var4d(outfid1,datashape4d,lonlength,latlength,altlength,timelength,&
721                 "fftau     ", "FFT ampl of vitu    ","m s-1     ",miss_val,&
722                  fftau )
723
724call write_var4d(outfid1,datashape4d,lonlength,latlength,altlength,timelength,&
725                 "ulf       ", "low freq part vitu  ","m s-1     ",miss_val,&
726                  ulf )
727
728call write_var4d(outfid1,datashape4d,lonlength,latlength,altlength,timelength,&
729                 "ubf       ", "band freq part vitu ","m s-1     ",miss_val,&
730                  ubf )
731
732call write_var4d(outfid1,datashape4d,lonlength,latlength,altlength,timelength,&
733                 "uhf       ", "high freq part vitu ","m s-1     ",miss_val,&
734                  uhf )
735endif !ok_out(1)
736
737if (ok_out(2)) then
738
739datashape4d(1)=lon_dimid2
740datashape4d(2)=lat_dimid2
741datashape4d(3)=alt_dimid2
742datashape4d(4)=time_dimid2
743datashape1d   =time_dimid2
744
745call write_var4d(outfid2,datashape4d,lonlength,latlength,altlength,timelength,&
746                 "fftav     ", "FFT ampl of vitv    ","m s-1     ",miss_val,&
747                  fftav )
748
749call write_var4d(outfid2,datashape4d,lonlength,latlength,altlength,timelength,&
750                 "vlf       ", "low freq part vitv  ","m s-1     ",miss_val,&
751                  vlf )
752
753call write_var4d(outfid2,datashape4d,lonlength,latlength,altlength,timelength,&
754                 "vbf       ", "band freq part vitv ","m s-1     ",miss_val,&
755                  vbf )
756
757call write_var4d(outfid2,datashape4d,lonlength,latlength,altlength,timelength,&
758                 "vhf       ", "high freq part vitv ","m s-1     ",miss_val,&
759                  vhf )
760endif !ok_out(2)
761
762if (ok_out(3)) then
763
764datashape4d(1)=lon_dimid3
765datashape4d(2)=lat_dimid3
766datashape4d(3)=alt_dimid3
767datashape4d(4)=time_dimid3
768datashape1d   =time_dimid3
769
770call write_var4d(outfid3,datashape4d,lonlength,latlength,altlength,timelength,&
771                 "fftaw     ", "FFT ampl of vitw    ","Pa s-1    ",miss_val,&
772                  fftaw )
773
774call write_var4d(outfid3,datashape4d,lonlength,latlength,altlength,timelength,&
775                 "wlf       ", "low freq part vitw  ","Pa s-1    ",miss_val,&
776                  wlf )
777
778call write_var4d(outfid3,datashape4d,lonlength,latlength,altlength,timelength,&
779                 "wbf       ", "band freq part vitw ","Pa s-1    ",miss_val,&
780                  wbf )
781
782 call write_var4d(outfid3,datashape4d,lonlength,latlength,altlength,timelength,&
783                 "whf       ", "high freq part vitw ","Pa s-1    ",miss_val,&
784                  whf )
785endif !ok_out(3)
786
787if (ok_out(4)) then
788
789datashape4d(1)=lon_dimid4
790datashape4d(2)=lat_dimid4
791datashape4d(3)=alt_dimid4
792datashape4d(4)=time_dimid4
793datashape1d   =time_dimid4
794
795call write_var4d(outfid4,datashape4d,lonlength,latlength,altlength,timelength,&
796                 "fftaT     ", "FFT ampl of temp    ","K         ",miss_val,&
797                  fftaT )
798
799call write_var4d(outfid4,datashape4d,lonlength,latlength,altlength,timelength,&
800                 "tlf       ", "low freq part temp  ","K         ",miss_val,&
801                  Tlf )
802
803call write_var4d(outfid4,datashape4d,lonlength,latlength,altlength,timelength,&
804                 "tbf       ", "band freq part temp ","K         ",miss_val,&
805                  Tbf )
806
807call write_var4d(outfid4,datashape4d,lonlength,latlength,altlength,timelength,&
808                 "thf       ", "high freq part temp ","K         ",miss_val,&
809                  Thf )
810endif !ok_out(4)
811
812!!!! Close output files
813if (ok_out(1)) then
814ierr=NF_CLOSE(outfid1)
815if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile1
816endif !ok_out(1)
817
818if (ok_out(2)) then
819ierr=NF_CLOSE(outfid2)
820if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile2
821endif !ok_out(2)
822
823if (ok_out(3)) then
824ierr=NF_CLOSE(outfid3)
825if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile3
826endif !ok_out(3)
827
828if (ok_out(4)) then
829ierr=NF_CLOSE(outfid4)
830if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile4
831endif !ok_out(4)
832
833
834end program
Note: See TracBrowser for help on using the repository browser.