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

Last change on this file since 3315 was 1454, checked in by slebonnois, 9 years ago

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