Changeset 1454 for trunk/LMDZ.TITAN/Tools/fft.F90

Timestamp:

Jun 16, 2015, 7:01:30 PM (9 years ago)

Author:

slebonnois

Message:

SL: corrections in Titan and Venus tools

File:

• trunk/LMDZ.TITAN/Tools/fft.F90 (modified) (10 diffs)

trunk/LMDZ.TITAN/Tools/fft.F90

@@ -63 +63 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -125 +124 @@
 logical :: lmdflag
 
+! Tuning parameters
+real :: fcoup1,fcoup2,width
+real :: fcoup1tmp,fcoup2tmp,widthtmp
+logical,dimension(4) :: ok_out
+character (len=1) :: ok_outtmp
+
 include "planet.h"
-include "filter.h"
 
 #include <fftw3.f> 
@@ -135 +139 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
 write(*,*) "You are working on the atmosphere of ",planet
+
+! initialisation
+!----------------
+
+! Par defaut
+
+! Define the filters
+! Low  cutting frequency, in Hz    : fcoup1
+fcoup1=2.5e-6
+! High cutting frequency, in Hz    : fcoup2
+fcoup2=6.5e-6
+! Half-width of the filters, in Hz : width
+width=4.e-7
+! Outputs (U,     V,      W,     T)
+ok_out=(/.true.,.true.,.false.,.true./)
+
+print*,"Low  cutting frequency, in Hz ?"
+print*,"between 1e-5 and 1e-7, 0 for default => 2.5e-6"
+read(*,*) fcoup1tmp
+if ((fcoup1tmp.lt.1e-5).and.(fcoup1tmp.gt.1e-7)) fcoup1=fcoup1tmp
+print*,"=",fcoup1
+print*,"High cutting frequency, in Hz ?"
+print*,"between 1e-5 and 1e-7, 0 for default => 6.5e-6"
+read(*,*) fcoup2tmp
+if ((fcoup2tmp.lt.1e-5).and.(fcoup2tmp.gt.1e-7)) fcoup2=fcoup2tmp
+print*,"=",fcoup2
+print*,"Half-width of the filters, in Hz ?"
+print*,"between 1e-6 and 1e-8, 0 for default => 4e-7)"
+read(*,*) widthtmp
+if ((widthtmp.lt.1e-6).and.(widthtmp.gt.1e-8)) width=widthtmp
+print*,"=",width
+!width = 3./time(timelength)
+
+! Outputs 
+print*,"Output of zonal wind ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(1)=.false.
+print*,"=",ok_out(1)
+print*,"Output of meridional wind ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(2)=.false.
+print*,"=",ok_out(2)
+print*,"Output of vertical wind ? (y or n, default is n)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."y") ok_out(3)=.true.
+print*,"=",ok_out(3)
+print*,"Output of temperature ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(4)=.false.
+print*,"=",ok_out(4)
 
 !===============================================================================
@@ -172 +227 @@
 ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
-
-allocate(latrad(latlength))
-latrad = lat*pi/180.
 
 allocate(plev(altlength))
@@ -212 +264 @@
 
 text="temp"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) then
     print*,"Error: Failed to get temperature ID"
@@ -236 +288 @@
 
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitu ID"
@@ -252 +304 @@
 
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitv ID"
@@ -268 +320 @@
 
 text="vitw"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitw ID"
@@ -288 +340 @@
 ! 2.2.1 FFT and filtering
 !===============================================================================
+
 ! allocations
 !-------------
@@ -405 +458 @@
 enddo
 
-! Define the filters
-
-! IN filter.h :
-! Low  cutting frequency, in Hz    : fcoup1
-! High cutting frequency, in Hz    : fcoup2
-! Half-width of the filters, in Hz : width
-
-!print*,"Low  cutting frequency, in Hz ?"
-!read(*,*) fcoup1
-!print*,"High cutting frequency, in Hz ?"
-!read(*,*) fcoup2
-!print*,"Half-width of the filters, in Hz ?"
-!read(*,*) width
-!width = 3./time(timelength)
-
 do itim=1,M_fft+1
   if (freq(itim).lt.(fcoup1-width)) then