Changeset 1426 for trunk/UTIL/SPECTRA

Timestamp:

May 11, 2015, 4:25:55 PM (10 years ago)

Author:

milmd

Message:

UTIL/SPECTRA. Correction of some bugs with spectra average. Now tested on Ciclad.

Location:

trunk/UTIL/SPECTRA

Files:

• readme (modified) (2 diffs)
• spectra_analysis.f90 (modified) (5 diffs)

trunk/UTIL/SPECTRA/readme

@@ -51 +51 @@
 FC=pgf90 or ifort
 
+for some machine, may load netcdf module:
+module load netcdf
+
+for ciclad machine, if using netcdf4:
+in makefile change line 7
+LDFLAGS=-L${netcdfpath}/lib -lnetcdf -L${spherepackpath}/lib -lspherepack
+and replace by
+LDFLAGS=-L${netcdfpath}/lib -lnetcdf -lnetcdff -L${spherepackpath}/lib -lspherepack 
+
 make spectra_analysis
 
@@ -56 +65 @@
 
 > compile test program
-make test_analysis
+make test_harmonic
 
 > create harmonic winds
-./test_analysis
+./test_harmonic
 
-> compute kinetic energy spectrum for one harmonic
+>>>>> some examples
+
+> compute kinetic energy spectrum for one harmonic test file
+
 ./spectra_analysis harmonic_64x48_lmdz_22.nc -alt none -time none -o harmonic_64x48_lmdz_22_spectra
 
+output file harmonic_64x48_lmdz_22_spectra contains:
+first column -> harmonic numbers 
+second column -> amplitude of harmonic component, here only component 22 is non zero.
+
+
+> compute kinetic energy spectrum with temporal and vertical mean
+
+./spectra_analysis diagfi.nc -t 5 -mt 10 -z 15 -mz 2 -o spectra
+
+compute kinetic energy spectrum for t=5,6,...,15 and z=15,16,17 and average all spectra obtained.
+
+
+> compute divergence and rotational part of kinetic energy spectrum
+
+./spectra_analysis diagfi.nc -t 5 -z 15 -o spectra -divrot
+
+in spectra file, will find the rotationnal and divergent part of the decomposition of the velocity on the vectorial spherical harmonic basis:
+#Spherical  diagfi.nc                                         
+#wavenumber t=  150 z=    1                                   
+#           spec_tot        spec_div        spec_rot          
+    0       0.000000E+00    0.000000E+00    0.000000E+00     
+    1       0.119347E+01    0.101938E+01    0.174095E+00     
+    2       0.530388E+00    0.250037E+00    0.280350E+00     
+    3       0.106973E+01    0.589848E+00    0.479880E+00     
+...
+
+
+>>>>> output interpretation
+
+One can plot velocity projection against spherical wavenumber with gnuplot for instance.
+In geostrophic turbulence, a n^-3 slope must appear for wavenumber n=1..100(?) and a n^-(5/3) slope for higher n.
+
+
+>>>>> reference
+J. N. Koshyk, 2001, The Horizontal Kinetic Energy Spectrum and Spectral Budget Simulated by a High-Resolution Troposphere–Stratosphere–Mesosphere GCM
+
+
+
+
+
+

trunk/UTIL/SPECTRA/spectra_analysis.f90

@@ -124 +124 @@
       i = i + 1
     elseif (arg(i) == '-h' .or. arg(i) == '--help') then
-      print*,'Usage\n spectra netcdfFile [option]\n [-h or --help]\t: brief help'
-      print*,'[-t int]\t: temporal indice (default: 1)\n [-z int]\t: vertical indice (default: 1)'
-      print*,'[-mt int]\t: extent of temporal mean (default: 0)\n [-mz int]\t: extent of vertical mean (default: 0)'
-      print*,'[-o str]\t: output spectra file name (default: spectra)'
-      print*,'[-u str]\t: name of zonal wind in netcdf file'
-      print*,'[-v str]\t: name of meridional wind in netcdf file'
-      print*,'[-lat str]\t: name of latitude field in netcdf file'
-      print*,'[-lon str]\t: name of longitude field in netcdf file'
-      print*,'[-alt str]\t: name of altitude field in netcdf file. ''none'' if no altitude axis.'
-      print*,'[-time str]\t: name of time field in netcdf file. ''none'' if no time axis.'
-      print*,'[-divrot]\t: total, divergence and vorticity spectra coefficient in output file'
+      print*,'Usage'
+      print*,'spectra_analysis netcdfFile [option]'
+      print*,'[-h or --help]    : brief help'
+      print*,'[-t int]  : temporal indice (default: 1)'
+      print*,'[-z int]  : vertical indice (default: 1)'
+      print*,'[-mt int] : extent of temporal mean (default: 0)'
+      print*,'[-mz int] : extent of vertical mean (default: 0)'
+      print*,'[-o str]  : output spectra file name (default: spectra)'
+      print*,'[-u str]  : name of zonal wind in netcdf file'
+      print*,'[-v str]  : name of meridional wind in netcdf file'
+      print*,'[-lat str]        : name of latitude field in netcdf file'
+      print*,'[-lon str]        : name of longitude field in netcdf file'
+      print*,'[-alt str]        : name of altitude field in netcdf file. ''none'' if no altitude axis.'
+      print*,'[-time str]       : name of time field in netcdf file. ''none'' if no time axis.'
+      print*,'[-divrot] : total, divergence and vorticity spectra coefficient in output file'
+      print*,''
+      print*,'Compute the kinetic energy spectrum of a wind field on a longitude-latitude grid.'
+      print*,'The grid must have a redundant point in longitude.'
+      print*,'Ideally the analysis works better on a symetric grid (2N points in longitude and N points in latitude).'
+      print*,'The output text file (called spectra by default) give the decomposition'
+      print*,'of the velocity on the vectorial spherical harmonic basis.'
       stop 'End help'
     else
@@ -475 +485 @@
 spectra_c_global(:,number_global) = spectra_c(:)
 
+!**********Some cleaning
+deallocate(br)
+deallocate(bi)
+deallocate(cr)
+deallocate(ci)
+deallocate(wvhaec)
+deallocate(dwork)
+deallocate(work)
+deallocate(norm)
+deallocate(norm_b)
+deallocate(norm_c)
+deallocate(spectra)
+deallocate(spectra_b)
+deallocate(spectra_c)
+deallocate(zonal_wind)
+deallocate(merid_wind)
+deallocate(zonal_wind_lmdz)
+deallocate(merid_wind_lmdz)
+
 !**********
 !End of loop over time and altitude indices
@@ -516 +545 @@
 end do
 end do
-write(10,*) ''
+write(10,*)
 write(10,'(a11,a1)',advance='no') '#wavenumber',' '
 do t = 1,n_indt
@@ -523 +552 @@
 end do
 end do
-write(10,*) ''
+write(10,*)
 write(10,'(a1,a11)',advance='no') '#',' '
 do t = 1,n_indt
@@ -534 +563 @@
 end do
 end do
-write(10,*)''
+write(10,*)
 !**********Writing file_spectra
 do j=1,nLat
-  write(10,'(i5,a7)',advance='no') j-1,' '
+  write(10,'(i5,a6)',advance='no') j-1,' '
   do t = 1,n_indt
   do z = 1,n_indz
   number_config = (t-1)*n_indz+z
   if (is_div_rot) then
-    write(10,'(e13.6E2,a4,e13.6E2,a4,e13.6E2,a6)',advance='no') spectra_config(j,number_config),' ',&
+    write(10,'(e13.6E2,a3,e13.6E2,a3,e13.6E2,a5)',advance='no') spectra_config(j,number_config),' ',&
                         spectra_b_config(j,number_config),' ',spectra_c_config(j,number_config),' '
   else
-    write(10,'(e13.6E2,a38)',advance='no') spectra_config(j,number_config),' '
+    write(10,'(e13.6E2,a37)',advance='no') spectra_config(j,number_config),' '
   end if
   end do
   end do
-  if (j /= nlat) write(10,*)''
+  if (j /= nlat) write(10,*)
 end do
 close(10) 
 
 print *, "***SUCCESS writing ",trim(file_spectra)
-deallocate(br)
-deallocate(bi)
-deallocate(cr)
-deallocate(ci)
-deallocate(wvhaec)
-deallocate(dwork)
-deallocate(work)
+
+!***********Some cleaning
+deallocate(spectra_config)
+deallocate(spectra_b_config)
+deallocate(spectra_c_config)
+deallocate(spectra_global)
+deallocate(spectra_b_global)
+deallocate(spectra_c_global)
+deallocate(tab_indt)
+deallocate(tab_indz)
+deallocate(tmp_tab_int)
+deallocate(arg)
 
 contains