source: trunk/LMDZ.VENUS/Tools/fft.F90 @ 937

Last change on this file since 937 was 888, checked in by slebonnois, 12 years ago

SL: small modifications to the tools, to Venus default .def files and to outputs (including forgotten modifications linked to the 1D); + bug corrections in phytitan
File size: 29.5 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!===============================================================================
210if (ok_out(4)) then
211allocate(temp(lonlength,latlength,altlength,timelength))
212
213text="temp"
214call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
215if (ierr1.ne.NF_NOERR) then
216  write(*,*) "  looking for t instead... "
217  text="t"
218  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
219  if (ierr1.ne.NF_NOERR) then
220    print*,"Error: Failed to get temperature ID"
221    ok_out(4)=.false.
222  endif
223endif
224if (ierr2.ne.NF_NOERR) then
225  print*,"Error: Failed reading temperature"
226  ok_out(4)=.false.
227endif
228endif !ok_out(4)
229
230!===============================================================================
231! 2.1.2 Winds
232!===============================================================================
233! zonal wind vitu (in m/s)
234if (ok_out(1)) then
235allocate(vitu(lonlength,latlength,altlength,timelength))
236
237text="vitu"
238call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
239if (ierr1.ne.NF_NOERR) then
240  print*,"Error: Failed to get vitu ID"
241  ok_out(1)=.false.
242endif
243if (ierr2.ne.NF_NOERR) then
244  print*,"Error: Failed reading zonal wind"
245  ok_out(1)=.false.
246endif
247endif !ok_out(1)
248
249! meridional wind vitv (in m/s)
250if (ok_out(2)) then
251allocate(vitv(lonlength,latlength,altlength,timelength))
252
253text="vitv"
254call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
255if (ierr1.ne.NF_NOERR) then
256  print*,"Error: Failed to get vitv ID"
257  ok_out(2)=.false.
258endif
259if (ierr2.ne.NF_NOERR) then
260  print*,"Error: Failed reading meridional wind"
261  ok_out(2)=.false.
262endif
263endif !ok_out(2)
264
265! vertical wind vitw (in Pa/s)
266if (ok_out(3)) then
267allocate(vitw(lonlength,latlength,altlength,timelength))
268
269text="vitw"
270call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
271if (ierr1.ne.NF_NOERR) then
272  print*,"Error: Failed to get vitw ID"
273  ok_out(3)=.false.
274endif
275if (ierr2.ne.NF_NOERR) then
276  print*,"Error: Failed reading vertical wind"
277  ok_out(3)=.false.
278endif
279endif !ok_out(3)
280
281!===============================================================================
282! 2.2 Computations
283!===============================================================================
284
285print*,"debut calcul"
286
287!===============================================================================
288! 2.2.1 FFT and filtering
289!===============================================================================
290! allocations
291!-------------
292if (ok_out(1)) then
293allocate(fftau(lonlength,latlength,altlength,timelength))
294allocate(uprim(lonlength,latlength,altlength,timelength))
295allocate(ulf(lonlength,latlength,altlength,timelength))
296allocate(ubf(lonlength,latlength,altlength,timelength))
297allocate(uhf(lonlength,latlength,altlength,timelength))
298endif !ok_out(1)
299if (ok_out(2)) then
300allocate(fftav(lonlength,latlength,altlength,timelength))
301allocate(vprim(lonlength,latlength,altlength,timelength))
302allocate(vlf(lonlength,latlength,altlength,timelength))
303allocate(vbf(lonlength,latlength,altlength,timelength))
304allocate(vhf(lonlength,latlength,altlength,timelength))
305endif !ok_out(2)
306if (ok_out(3)) then
307allocate(fftaw(lonlength,latlength,altlength,timelength))
308allocate(wprim(lonlength,latlength,altlength,timelength))
309allocate(wlf(lonlength,latlength,altlength,timelength))
310allocate(wbf(lonlength,latlength,altlength,timelength))
311allocate(whf(lonlength,latlength,altlength,timelength))
312endif !ok_out(3)
313if (ok_out(4)) then
314allocate(fftaT(lonlength,latlength,altlength,timelength))
315allocate(Tprim(lonlength,latlength,altlength,timelength))
316allocate(Tlf(lonlength,latlength,altlength,timelength))
317allocate(Tbf(lonlength,latlength,altlength,timelength))
318allocate(Thf(lonlength,latlength,altlength,timelength))
319endif !ok_out(4)
320
321! lon,lat,alt
322if (ok_out(1)) allocate(umean(lonlength,latlength,altlength))
323if (ok_out(2)) allocate(vmean(lonlength,latlength,altlength))
324if (ok_out(3)) allocate(wmean(lonlength,latlength,altlength))
325if (ok_out(4)) allocate(Tmean(lonlength,latlength,altlength))
326
327! time / frequencies
328allocate(freq(timelength))
329allocate(wndow(timelength))
330allocate(var(timelength))
331allocate(fltvar(timelength))
332M_fft = timelength/2
333allocate(fftvar(M_fft+1))
334allocate(fltfft(M_fft+1))
335allocate(filtrelf(M_fft+1))
336allocate(filtrebf(M_fft+1))
337allocate(filtrehf(M_fft+1))
338
339! intermediates
340!-----------------
341
342if (ok_out(1)) call moytim(lonlength,latlength,altlength,timelength,miss_val,vitu,umean)
343if (ok_out(2)) call moytim(lonlength,latlength,altlength,timelength,miss_val,vitv,vmean)
344if (ok_out(3)) call moytim(lonlength,latlength,altlength,timelength,miss_val,vitw,wmean)
345if (ok_out(4)) call moytim(lonlength,latlength,altlength,timelength,miss_val,temp,Tmean)
346
347do ilon=1,lonlength
348 do ilat=1,latlength
349  do ilev=1,altlength
350   do itim=1,timelength
351if (ok_out(1)) then
352    if ((vitu(ilon,ilat,ilev,itim).ne.miss_val).and. &
353        (umean(ilon,ilat,ilev)    .ne.miss_val)) then
354  uprim(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim)-umean(ilon,ilat,ilev)
355    else
356  uprim(ilon,ilat,ilev,itim) = miss_val
357    endif
358endif !ok_out(1)
359if (ok_out(2)) then
360    if ((vitv(ilon,ilat,ilev,itim).ne.miss_val).and. &
361        (vmean(ilon,ilat,ilev)    .ne.miss_val)) then
362  vprim(ilon,ilat,ilev,itim) = vitv(ilon,ilat,ilev,itim)-vmean(ilon,ilat,ilev)
363    else
364  vprim(ilon,ilat,ilev,itim) = miss_val
365    endif
366endif !ok_out(2)
367if (ok_out(3)) then
368    if ((vitw(ilon,ilat,ilev,itim).ne.miss_val).and. &
369        (wmean(ilon,ilat,ilev)    .ne.miss_val)) then
370  wprim(ilon,ilat,ilev,itim) = vitw(ilon,ilat,ilev,itim)-wmean(ilon,ilat,ilev)
371    else
372  wprim(ilon,ilat,ilev,itim) = miss_val
373    endif
374endif !ok_out(3)
375if (ok_out(4)) then
376    if ((temp(ilon,ilat,ilev,itim).ne.miss_val).and. &
377        (Tmean(ilon,ilat,ilev)    .ne.miss_val)) then
378  Tprim(ilon,ilat,ilev,itim) = temp(ilon,ilat,ilev,itim)-Tmean(ilon,ilat,ilev)
379    else
380  Tprim(ilon,ilat,ilev,itim) = miss_val
381    endif
382endif !ok_out(4)
383   enddo
384  enddo
385 enddo
386enddo ! lonlength
387
388! fft intermediates
389!-------------
390
391! Define the frequencies
392do itim=1,M_fft+1
393  freq(itim) = (itim-1)/(timelength*(time(2)-time(1)))
394enddo
395do itim=M_fft+2,timelength
396  freq(itim) = 0.
397enddo
398
399! Define the window (triangle)
400do itim=1,timelength
401! N window:
402!  wndow(itim)= 1.
403! triangulaire de moyenne = 1
404  wndow(itim)= 2.*(1. - abs(real(itim-0.5-M_fft)/real(M_fft)))
405enddo
406
407! Define the filters
408
409! IN filter.h :
410! Low  cutting frequency, in Hz    : fcoup1
411! High cutting frequency, in Hz    : fcoup2
412! Half-width of the filters, in Hz : width
413
414!print*,"Low  cutting frequency, in Hz ?"
415!read(*,*) fcoup1
416!print*,"High cutting frequency, in Hz ?"
417!read(*,*) fcoup2
418!print*,"Half-width of the filters, in Hz ?"
419!read(*,*) width
420!width = 3./time(timelength)
421
422do itim=1,M_fft+1
423  if (freq(itim).lt.(fcoup1-width)) then
424     filtrelf(itim) = 1.
425  elseif (freq(itim).gt.(fcoup1+width)) then
426     filtrelf(itim) = 0.
427  else
428     filtrelf(itim) = (1.+sin(pi*(fcoup1-freq(itim))/(2.*width)))/2.
429  endif
430  if (freq(itim).lt.(fcoup2-width)) then
431     filtrehf(itim) = 0.
432  elseif (freq(itim).gt.(fcoup2+width)) then
433     filtrehf(itim) = 1.
434  else
435     filtrehf(itim) = (1.-sin(pi*(fcoup2-freq(itim))/(2.*width)))/2.
436  endif
437  filtrebf(itim) = (1.-filtrelf(itim))*(1.-filtrehf(itim))
438enddo
439
440
441! fft and filtering
442!-------------
443
444!---FFTW routines
445call dfftw_plan_dft_r2c_1d(planf,timelength,var,fftvar,FFTW_MEASURE)
446call dfftw_plan_dft_c2r_1d(planb,timelength,fltfft,fltvar,FFTW_MEASURE)
447!---
448
449do ilon=1,lonlength
450 do ilat=1,latlength
451  do ilev=1,altlength
452
453! For zonal wind field
454if (ok_out(1)) then
455
456   flagfft=.true.
457   do itim=1,timelength
458     if (uprim(ilon,ilat,ilev,itim).eq.miss_val) flagfft=.false.
459   enddo
460
461   if (flagfft) then
462
463! 1/ windowing to improve spectral analysis
464      var(:)=uprim(ilon,ilat,ilev,:)*wndow(:)
465! 2/ FFT computation
466!---FFTW routines
467      call dfftw_execute_dft_r2c(planf,var,fftvar)
468!---
469! 3/ Amplitude of the FFT, for spectral analysis
470      fftau(ilon,ilat,ilev,1)=abs(fftvar(1))/M_fft
471      do itim=2,M_fft
472       fftau(ilon,ilat,ilev,itim) = abs(fftvar(itim))/M_fft
473      enddo
474      fftau(ilon,ilat,ilev,M_fft+1)=abs(fftvar(M_fft+1))/M_fft
475      do itim=M_fft+2,timelength
476       fftau(ilon,ilat,ilev,itim) = 0.
477      enddo
478
479! 4/ filtering the FFT in three regions
480! filtering + normalisation (low freq)
481      fltfft(:) = fftvar(:)*filtrelf(:)/timelength
482! 5/ backward FFT for each region
483!---FFTW routines
484      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
485!---
486! 6/ reverse the windowing
487      ulf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
488
489! filtering + normalisation (band freq)
490      fltfft(:) = fftvar(:)*filtrebf(:)/timelength
491!---FFTW routines
492      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
493!---
494      ubf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
495
496! filtering + normalisation (high freq)
497      fltfft(:) = fftvar(:)*filtrehf(:)/timelength
498!---FFTW routines
499      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
500!---
501      uhf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
502
503   else
504     fftau(ilon,ilat,ilev,itim) = miss_val
505       ulf(ilon,ilat,ilev,itim) = miss_val
506       ubf(ilon,ilat,ilev,itim) = miss_val
507       uhf(ilon,ilat,ilev,itim) = miss_val
508   endif ! flagfft
509
510endif !ok_out(1)
511
512! For meridional wind wind field
513if (ok_out(2)) then
514
515   flagfft=.true.
516   do itim=1,timelength
517     if (vprim(ilon,ilat,ilev,itim).eq.miss_val) flagfft=.false.
518   enddo
519
520   if (flagfft) then
521
522! 1/ windowing to improve spectral analysis
523      var(:)=vprim(ilon,ilat,ilev,:)*wndow(:)
524! 2/ FFT computation
525!---FFTW routines
526      call dfftw_execute_dft_r2c(planf,var,fftvar)
527!---
528! 3/ Amplitude of the FFT, for spectral analysis
529      fftav(ilon,ilat,ilev,1)=abs(fftvar(1))/M_fft
530      do itim=2,M_fft
531       fftav(ilon,ilat,ilev,itim) = abs(fftvar(itim))/M_fft
532      enddo
533      fftav(ilon,ilat,ilev,M_fft+1)=abs(fftvar(M_fft+1))/M_fft
534      do itim=M_fft+2,timelength
535       fftav(ilon,ilat,ilev,itim) = 0.
536      enddo
537
538! 4/ filtering the FFT in three regions
539! filtering + normalisation (low freq)
540      fltfft(:) = fftvar(:)*filtrelf(:)/timelength
541! 5/ backward FFT for each region
542!---FFTW routines
543      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
544!---
545! 6/ reverse the windowing
546      vlf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
547
548! filtering + normalisation (band freq)
549      fltfft(:) = fftvar(:)*filtrebf(:)/timelength
550!---FFTW routines
551      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
552!---
553      vbf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
554
555! filtering + normalisation (high freq)
556      fltfft(:) = fftvar(:)*filtrehf(:)/timelength
557!---FFTW routines
558      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
559!---
560      vhf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
561
562   else
563     fftav(ilon,ilat,ilev,itim) = miss_val
564       vlf(ilon,ilat,ilev,itim) = miss_val
565       vbf(ilon,ilat,ilev,itim) = miss_val
566       vhf(ilon,ilat,ilev,itim) = miss_val
567   endif ! flagfft
568
569endif !ok_out(2)
570
571! For vertical wind wind field
572if (ok_out(3)) then
573
574   flagfft=.true.
575   do itim=1,timelength
576     if (wprim(ilon,ilat,ilev,itim).eq.miss_val) flagfft=.false.
577   enddo
578
579   if (flagfft) then
580
581! 1/ windowing to improve spectral analysis
582      var(:)=wprim(ilon,ilat,ilev,:)*wndow(:)
583! 2/ FFT computation
584!---FFTW routines
585      call dfftw_execute_dft_r2c(planf,var,fftvar)
586!---
587! 3/ Amplitude of the FFT, for spectral analysis
588      fftaw(ilon,ilat,ilev,1)=abs(fftvar(1))/M_fft
589      do itim=2,M_fft
590       fftaw(ilon,ilat,ilev,itim) = abs(fftvar(itim))/M_fft
591      enddo
592      fftaw(ilon,ilat,ilev,M_fft+1)=abs(fftvar(M_fft+1))/M_fft
593      do itim=M_fft+2,timelength
594       fftaw(ilon,ilat,ilev,itim) = 0.
595      enddo
596
597! 4/ filtering the FFT in three regions
598! filtering + normalisation (low freq)
599      fltfft(:) = fftvar(:)*filtrelf(:)/timelength
600! 5/ backward FFT for each region
601!---FFTW routines
602      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
603!---
604! 6/ reverse the windowing
605      wlf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
606
607! filtering + normalisation (band freq)
608      fltfft(:) = fftvar(:)*filtrebf(:)/timelength
609!---FFTW routines
610      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
611!---
612      wbf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
613
614! filtering + normalisation (high freq)
615      fltfft(:) = fftvar(:)*filtrehf(:)/timelength
616!---FFTW routines
617      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
618!---
619      whf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
620
621   else
622     fftaw(ilon,ilat,ilev,itim) = miss_val
623       wlf(ilon,ilat,ilev,itim) = miss_val
624       wbf(ilon,ilat,ilev,itim) = miss_val
625       whf(ilon,ilat,ilev,itim) = miss_val
626   endif ! flagfft
627
628endif !ok_out(3)
629
630! For temperature field
631if (ok_out(4)) then
632
633   flagfft=.true.
634   do itim=1,timelength
635     if (Tprim(ilon,ilat,ilev,itim).eq.miss_val) flagfft=.false.
636   enddo
637
638   if (flagfft) then
639
640! 1/ windowing to improve spectral analysis
641      var(:)=Tprim(ilon,ilat,ilev,:)*wndow(:)
642! 2/ FFT computation
643!---FFTW routines
644      call dfftw_execute_dft_r2c(planf,var,fftvar)
645!---
646! 3/ Amplitude of the FFT, for spectral analysis
647      fftaT(ilon,ilat,ilev,1)=abs(fftvar(1))/M_fft
648      do itim=2,M_fft
649       fftaT(ilon,ilat,ilev,itim) = abs(fftvar(itim))/M_fft
650      enddo
651      fftaT(ilon,ilat,ilev,M_fft+1)=abs(fftvar(M_fft+1))/M_fft
652      do itim=M_fft+2,timelength
653       fftaT(ilon,ilat,ilev,itim) = 0.
654      enddo
655
656! 4/ filtering the FFT in three regions
657! filtering + normalisation (low freq)
658      fltfft(:) = fftvar(:)*filtrelf(:)/timelength
659! 5/ backward FFT for each region
660!---FFTW routines
661      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
662!---
663! 6/ reverse the windowing
664      Tlf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
665
666! filtering + normalisation (band freq)
667      fltfft(:) = fftvar(:)*filtrebf(:)/timelength
668!---FFTW routines
669      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
670!---
671      Tbf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
672
673! filtering + normalisation (high freq)
674      fltfft(:) = fftvar(:)*filtrehf(:)/timelength
675!---FFTW routines
676      call dfftw_execute_dft_c2r(planb,fltfft,fltvar)
677!---
678      Thf(ilon,ilat,ilev,:) = fltvar(:)/wndow(:)
679
680   else
681     fftaT(ilon,ilat,ilev,itim) = miss_val
682       Tlf(ilon,ilat,ilev,itim) = miss_val
683       Tbf(ilon,ilat,ilev,itim) = miss_val
684       Thf(ilon,ilat,ilev,itim) = miss_val
685   endif ! flagfft
686
687endif !ok_out(4)
688
689  enddo
690 enddo
691enddo ! lonlength
692
693!---FFTW routines
694call dfftw_destroy_plan(planf)
695call dfftw_destroy_plan(planb)
696!---
697
698print*,"End of computations"
699
700!===============================================================================
701! 3. Create output files
702!===============================================================================
703
704! Create output files
705if (ok_out(1)) then
706ierr=NF_CREATE(outfile1,NF_CLOBBER,outfid1)
707if (ierr.ne.NF_NOERR) then
708  write(*,*)"Error: could not create file ",outfile1
709  stop
710endif
711endif !ok_out(1)
712
713if (ok_out(2)) then
714ierr=NF_CREATE(outfile2,NF_CLOBBER,outfid2)
715if (ierr.ne.NF_NOERR) then
716  write(*,*)"Error: could not create file ",outfile2
717  stop
718endif
719endif !ok_out(2)
720
721if (ok_out(3)) then
722ierr=NF_CREATE(outfile3,NF_CLOBBER,outfid3)
723if (ierr.ne.NF_NOERR) then
724  write(*,*)"Error: could not create file ",outfile3
725  stop
726endif
727endif !ok_out(3)
728
729if (ok_out(4)) then
730ierr=NF_CREATE(outfile4,NF_CLOBBER,outfid4)
731if (ierr.ne.NF_NOERR) then
732  write(*,*)"Error: could not create file ",outfile4
733  stop
734endif
735endif !ok_out(4)
736
737!===============================================================================
738! 3.1. Define and write dimensions
739!===============================================================================
740
741if (ok_out(1)) &
742call write_dim(outfid1,lonlength,latlength,altlength,timelength, &
743    lon,lat,plev,time,lon_dimid1,lat_dimid1,alt_dimid1,time_dimid1)
744if (ok_out(2)) &
745call write_dim(outfid2,lonlength,latlength,altlength,timelength, &
746    lon,lat,plev,time,lon_dimid2,lat_dimid2,alt_dimid2,time_dimid2)
747if (ok_out(3)) &
748call write_dim(outfid3,lonlength,latlength,altlength,timelength, &
749    lon,lat,plev,time,lon_dimid3,lat_dimid3,alt_dimid3,time_dimid3)
750if (ok_out(4)) &
751call write_dim(outfid4,lonlength,latlength,altlength,timelength, &
752    lon,lat,plev,time,lon_dimid4,lat_dimid4,alt_dimid4,time_dimid4)
753
754!===============================================================================
755! 3.2. Define and write variables
756!===============================================================================
757
758if (ok_out(1)) then
759
760datashape4d(1)=lon_dimid1
761datashape4d(2)=lat_dimid1
762datashape4d(3)=alt_dimid1
763datashape4d(4)=time_dimid1
764datashape1d   =time_dimid1
765
766call write_var1d(outfid1,datashape1d,timelength,&
767                "freq      ", "FFT frequencies     ","s-1       ",miss_val,&
768                 freq )
769
770call write_var4d(outfid1,datashape4d,lonlength,latlength,altlength,timelength,&
771                 "fftau     ", "FFT ampl of vitu    ","m s-1     ",miss_val,&
772                  fftau )
773
774call write_var4d(outfid1,datashape4d,lonlength,latlength,altlength,timelength,&
775                 "ulf       ", "low freq part vitu  ","m s-1     ",miss_val,&
776                  ulf )
777
778call write_var4d(outfid1,datashape4d,lonlength,latlength,altlength,timelength,&
779                 "ubf       ", "band freq part vitu ","m s-1     ",miss_val,&
780                  ubf )
781
782call write_var4d(outfid1,datashape4d,lonlength,latlength,altlength,timelength,&
783                 "uhf       ", "high freq part vitu ","m s-1     ",miss_val,&
784                  uhf )
785endif !ok_out(1)
786
787if (ok_out(2)) then
788
789datashape4d(1)=lon_dimid2
790datashape4d(2)=lat_dimid2
791datashape4d(3)=alt_dimid2
792datashape4d(4)=time_dimid2
793datashape1d   =time_dimid2
794
795call write_var1d(outfid2,datashape1d,timelength,&
796                "freq      ", "FFT frequencies     ","s-1       ",miss_val,&
797                 freq )
798
799call write_var4d(outfid2,datashape4d,lonlength,latlength,altlength,timelength,&
800                 "fftav     ", "FFT ampl of vitv    ","m s-1     ",miss_val,&
801                  fftav )
802
803call write_var4d(outfid2,datashape4d,lonlength,latlength,altlength,timelength,&
804                 "vlf       ", "low freq part vitv  ","m s-1     ",miss_val,&
805                  vlf )
806
807call write_var4d(outfid2,datashape4d,lonlength,latlength,altlength,timelength,&
808                 "vbf       ", "band freq part vitv ","m s-1     ",miss_val,&
809                  vbf )
810
811call write_var4d(outfid2,datashape4d,lonlength,latlength,altlength,timelength,&
812                 "vhf       ", "high freq part vitv ","m s-1     ",miss_val,&
813                  vhf )
814endif !ok_out(2)
815
816if (ok_out(3)) then
817
818datashape4d(1)=lon_dimid3
819datashape4d(2)=lat_dimid3
820datashape4d(3)=alt_dimid3
821datashape4d(4)=time_dimid3
822datashape1d   =time_dimid3
823
824call write_var1d(outfid3,datashape1d,timelength,&
825                "freq      ", "FFT frequencies     ","s-1       ",miss_val,&
826                 freq )
827
828call write_var4d(outfid3,datashape4d,lonlength,latlength,altlength,timelength,&
829                 "fftaw     ", "FFT ampl of vitw    ","Pa s-1    ",miss_val,&
830                  fftaw )
831
832call write_var4d(outfid3,datashape4d,lonlength,latlength,altlength,timelength,&
833                 "wlf       ", "low freq part vitw  ","Pa s-1    ",miss_val,&
834                  wlf )
835
836call write_var4d(outfid3,datashape4d,lonlength,latlength,altlength,timelength,&
837                 "wbf       ", "band freq part vitw ","Pa s-1    ",miss_val,&
838                  wbf )
839
840 call write_var4d(outfid3,datashape4d,lonlength,latlength,altlength,timelength,&
841                 "whf       ", "high freq part vitw ","Pa s-1    ",miss_val,&
842                  whf )
843endif !ok_out(3)
844
845if (ok_out(4)) then
846
847datashape4d(1)=lon_dimid4
848datashape4d(2)=lat_dimid4
849datashape4d(3)=alt_dimid4
850datashape4d(4)=time_dimid4
851datashape1d   =time_dimid4
852
853call write_var1d(outfid4,datashape1d,timelength,&
854                "freq      ", "FFT frequencies     ","s-1       ",miss_val,&
855                 freq )
856
857call write_var4d(outfid4,datashape4d,lonlength,latlength,altlength,timelength,&
858                 "fftaT     ", "FFT ampl of temp    ","K         ",miss_val,&
859                  fftaT )
860
861call write_var4d(outfid4,datashape4d,lonlength,latlength,altlength,timelength,&
862                 "tlf       ", "low freq part temp  ","K         ",miss_val,&
863                  Tlf )
864
865call write_var4d(outfid4,datashape4d,lonlength,latlength,altlength,timelength,&
866                 "tbf       ", "band freq part temp ","K         ",miss_val,&
867                  Tbf )
868
869call write_var4d(outfid4,datashape4d,lonlength,latlength,altlength,timelength,&
870                 "thf       ", "high freq part temp ","K         ",miss_val,&
871                  Thf )
872endif !ok_out(4)
873
874!!!! Close output files
875if (ok_out(1)) then
876ierr=NF_CLOSE(outfid1)
877if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile1
878endif !ok_out(1)
879
880if (ok_out(2)) then
881ierr=NF_CLOSE(outfid2)
882if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile2
883endif !ok_out(2)
884
885if (ok_out(3)) then
886ierr=NF_CLOSE(outfid3)
887if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile3
888endif !ok_out(3)
889
890if (ok_out(4)) then
891ierr=NF_CLOSE(outfid4)
892if (ierr.ne.NF_NOERR) write(*,*) 'Error, failed to close output file ',outfile4
893endif !ok_out(4)
894
895
896end program
Note: See TracBrowser for help on using the repository browser.