[cig-commits] [commit] Include_averaging_programs, master: Add time averaging programs (7ca97d3)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 4 20:00:21 PDT 2014
Repository : https://github.com/geodynamics/calypso
On branches: Include_averaging_programs,master
Link : https://github.com/geodynamics/calypso/compare/d3671756c9a8dc30cd0c89470f2e2f0599bbe718...1d20f55fd74b30f7715f875492f6a7b62f3b6274
>---------------------------------------------------------------
commit 7ca97d37e5511892724978f5d55f9f4d16aa968b
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Wed Apr 23 17:48:28 2014 -0700
Add time averaging programs
>---------------------------------------------------------------
7ca97d37e5511892724978f5d55f9f4d16aa968b
src/programs/Makefile | 2 +-
src/programs/data_utilities/Makefile | 3 +-
.../data_utilities/TIME_HISTORIES/Makefile | 67 ++++
.../TIME_HISTORIES/cal_tave_sph_ene_spectr.f90 | 248 ++++++++++++
.../TIME_HISTORIES/m_sph_ene_spectra.f90 | 442 +++++++++++++++++++++
.../TIME_HISTORIES/m_tave_sph_ene_spectr.f90 | 280 +++++++++++++
.../TIME_HISTORIES/t_average_sph_ene_spec.f90 | 123 ++++++
.../TIME_HISTORIES/tave_picked_sph_spec_data.f90 | 171 ++++++++
8 files changed, 1334 insertions(+), 2 deletions(-)
diff --git a/src/programs/Makefile b/src/programs/Makefile
index c09c2bd..0283bf2 100644
--- a/src/programs/Makefile
+++ b/src/programs/Makefile
@@ -52,7 +52,7 @@ target_task:
@echo 'parallels: sph_mhd' \
>> $(MAKENAME)
@echo '' >> $(MAKENAME)
- @echo 'utilities: data_utils mesh_utils' >> $(MAKENAME)
+ @echo 'utilities: data_utils mesh_utils sph_data_util' >> $(MAKENAME)
@echo '' >> $(MAKENAME)
@for dir in $(SUBDIRS); do \
( cd $${dir}; \
diff --git a/src/programs/data_utilities/Makefile b/src/programs/data_utilities/Makefile
index be9c44b..6be2ac8 100644
--- a/src/programs/data_utilities/Makefile
+++ b/src/programs/data_utilities/Makefile
@@ -5,7 +5,8 @@
SUBDIRS = \
MERGE \
INITIAL_FIELD \
-SNAPSHOT_MHD
+SNAPSHOT_MHD \
+TIME_HISTORIES
#
diff --git a/src/programs/data_utilities/TIME_HISTORIES/Makefile b/src/programs/data_utilities/TIME_HISTORIES/Makefile
new file mode 100644
index 0000000..259d683
--- /dev/null
+++ b/src/programs/data_utilities/TIME_HISTORIES/Makefile
@@ -0,0 +1,67 @@
+#
+#
+#
+
+SPH_UTILS_DIR = $$(DATA_UTILS_DIR)/TIME_HISTORIES
+
+TARGET_TAVE_SPH_ENE_SPEC = t_ave_sph_mean_square
+TARGET_T_AVE_PICK_SPH_SPEC = t_ave_picked_sph_coefs
+
+SOURCES = $(shell ls *.f90)
+
+MOD_TAVE_SPH_ENE_SPEC = \
+t_average_sph_ene_spec.o \
+m_sph_ene_spectra.o \
+m_tave_sph_ene_spectr.o \
+cal_tave_sph_ene_spectr.o
+
+MOD_T_AVE_PICK_SPH = \
+tave_picked_sph_spec_data.o
+
+#
+# ------------------------------------------------------------------
+#
+
+dir_list:
+ @echo 'SPH_UTILS_DIR = $(SPH_UTILS_DIR)' >> $(MAKENAME)
+
+target_list:
+ @echo 'TARGET_TAVE_SPH_ENE_SPEC = $$(BUILDDIR)/$(TARGET_TAVE_SPH_ENE_SPEC)' \
+ >> $(MAKENAME)
+ @echo 'TARGET_T_AVE_PICK_SPH_SPEC = $$(BUILDDIR)/$(TARGET_T_AVE_PICK_SPH_SPEC)' \
+ >> $(MAKENAME)
+ @echo >> $(MAKENAME)
+
+target_task:
+ @echo sph_data_util: \
+ '$$(TARGET_TAVE_SPH_ENE_SPEC) $$(TARGET_T_AVE_PICK_SPH_SPEC)' \
+ '$$(TARGET_MAX_SPH_ENE) ' >> $(MAKENAME)
+ @echo >> $(MAKENAME)
+ @echo '$$(TARGET_TAVE_SPH_ENE_SPEC)': '$$(MOD_TAVE_SPH_ENE_SPEC)' \
+ '$$(LIB_FILES_SPH_MHD)' \
+ >> $(MAKENAME)
+ @echo ' ''$$(F90)' '$$(F90FLAGS)' -o '$$(TARGET_TAVE_SPH_ENE_SPEC)' \
+ '$$(MOD_TAVE_SPH_ENE_SPEC)' \
+ '-L. $$(LIBS_SPH_MHD)' \
+ '$$(F90LIBS)' >> $(MAKENAME)
+ @echo '' >> $(MAKENAME)
+ @echo '$$(TARGET_T_AVE_PICK_SPH_SPEC)': '$$(MOD_T_AVE_PICK_SPH)' \
+ '$$(LIB_FILES_SPH_MHD)' \
+ >> $(MAKENAME)
+ @echo ' ''$$(F90)' '$$(F90FLAGS)' -o '$$(TARGET_T_AVE_PICK_SPH_SPEC)' \
+ '$$(MOD_T_AVE_PICK_SPH)' \
+ '-L. $$(LIBS_SPH_MHD)' \
+ '$$(F90LIBS)' >> $(MAKENAME)
+ @echo '' >> $(MAKENAME)
+
+
+lib_name:
+
+mod_list:
+ @echo MOD_TAVE_SPH_ENE_SPEC= \\ >> $(MAKENAME)
+ @echo $(MOD_TAVE_SPH_ENE_SPEC) >> $(MAKENAME)
+ @echo MOD_T_AVE_PICK_SPH= \\ >> $(MAKENAME)
+ @echo $(MOD_T_AVE_PICK_SPH) >> $(MAKENAME)
+
+module:
+ @$(MAKE_MOD_DEP) '$(MAKENAME)' '$$(SPH_UTILS_DIR)' $(SOURCES)
diff --git a/src/programs/data_utilities/TIME_HISTORIES/cal_tave_sph_ene_spectr.f90 b/src/programs/data_utilities/TIME_HISTORIES/cal_tave_sph_ene_spectr.f90
new file mode 100644
index 0000000..85f207e
--- /dev/null
+++ b/src/programs/data_utilities/TIME_HISTORIES/cal_tave_sph_ene_spectr.f90
@@ -0,0 +1,248 @@
+!>@file cal_tave_sph_ene_spectr.f90
+!! module cal_tave_sph_ene_spectr
+!!
+!! @author H. Matsui
+!! @date Programmed in May, 2008
+!!
+!
+!> @brief Evaluate time average spherical harmonics spectrum data
+!!
+!!@verbatim
+!! subroutine sum_average_ene_sph
+!! subroutine sum_deviation_ene_sph
+!! subroutine divide_average_ene_sph
+!! subroutine divide_deviation_ene_sph
+!!
+!! subroutine output_tave_ene_sph_data
+!! subroutine output_tsigma_ene_sph_data
+!!
+!! subroutine count_degree_on_layer_data
+!! subroutine count_degree_one_layer_data
+!! subroutine count_degree_on_volume_data
+!!@endverbatim
+!!
+!!@n @param istep time step number
+!!@n @param icou counter for snapshots
+!!@n @param ierr error flag
+!
+ module cal_tave_sph_ene_spectr
+!
+ use m_precision
+ use m_constants
+ use m_sph_ene_spectra
+ use m_tave_sph_ene_spectr
+!
+ implicit none
+!
+! --------------------------------------------------------------------
+!
+ contains
+!
+! --------------------------------------------------------------------
+!
+ subroutine sum_average_ene_sph
+!
+ integer(kind = kint) :: kr, nd, lth
+ real(kind= kreal) :: tmp, tmp_l, tmp_m, tmp_lm
+!
+!
+!$omp parallel private(kr,lth,nd,tmp,tmp_l,tmp_m,tmp_lm)
+ do kr = 1, nri_sph
+!$omp do
+ do nd = 1, ncomp_sph_spec
+ tmp = spectr_t(nd,kr)
+ ave_spec_t(nd,kr) = ave_spec_t(nd,kr) &
+ & + half * (tmp + spectr_pre_t(nd,kr)) &
+ & * (time_sph + pre_time)
+ spectr_pre_t(nd,kr) = tmp
+ end do
+!$omp end do nowait
+!
+ do lth = 0, ltr_sph
+!$omp do
+ do nd = 1, ncomp_sph_spec
+ tmp_l = spectr_l(nd,lth,kr)
+ tmp_m = spectr_m(nd,lth,kr)
+ tmp_lm = spectr_lm(nd,lth,kr)
+!
+ ave_spec_l(nd,lth,kr) = ave_spec_l(nd,lth,kr) &
+ & + half * (tmp_l + spectr_pre_l(nd,lth,kr)) &
+ & * (time_sph + pre_time)
+ ave_spec_m(nd,lth,kr) = ave_spec_m(nd,lth,kr) &
+ & + half * (tmp_m + spectr_pre_l(nd,lth,kr)) &
+ & * (time_sph + pre_time)
+ ave_spec_lm(nd,lth,kr) = ave_spec_lm(nd,lth,kr) &
+ & + half * (tmp_lm + spectr_pre_l(nd,lth,kr)) &
+ & * (time_sph + pre_time)
+!
+ spectr_pre_l(nd,lth,kr) = tmp_l
+ spectr_pre_m(nd,lth,kr) = tmp_m
+ spectr_pre_lm(nd,lth,kr) = tmp_lm
+ end do
+!$omp end do nowait
+ end do
+ end do
+!$omp end parallel
+!
+ pre_time = time_sph
+!
+ end subroutine sum_average_ene_sph
+!
+! --------------------------------------------------------------------
+!
+ subroutine sum_deviation_ene_sph
+!
+ integer(kind = kint) :: kr, nd, lth
+ real(kind= kreal) :: tmp, tmp_l, tmp_m, tmp_lm
+!
+!
+!$omp parallel private(kr,lth,nd,tmp,tmp_l,tmp_m,tmp_lm)
+ do kr = 1, nri_sph
+!$omp do
+ do nd = 1, ncomp_sph_spec
+ tmp = (spectr_t(nd,kr) - ave_spec_t(nd,kr))**2
+ sigma_spec_t(nd,kr) = sigma_spec_t(nd,kr) &
+ & + half * (tmp + spectr_pre_t(nd,kr)) &
+ & * (time_sph + pre_time)
+ spectr_pre_t(nd,kr) = tmp
+ end do
+!$omp end do nowait
+!
+ do lth = 0, ltr_sph
+!$omp do
+ do nd = 1, ncomp_sph_spec
+ tmp_l = (spectr_l(nd,lth,kr) - ave_spec_l(nd,lth,kr))**2
+ tmp_m = (spectr_m(nd,lth,kr) - ave_spec_m(nd,lth,kr))**2
+ tmp_lm = (spectr_lm(nd,lth,kr) - ave_spec_lm(nd,lth,kr))**2
+!
+ sigma_spec_l(nd,lth,kr) = sigma_spec_l(nd,lth,kr) &
+ & + half * (tmp_l + spectr_pre_l(nd,lth,kr)) &
+ & * (time_sph + pre_time)
+ sigma_spec_m(nd,lth,kr) = sigma_spec_m(nd,lth,kr) &
+ & + half * (tmp_m + spectr_pre_l(nd,lth,kr)) &
+ & * (time_sph + pre_time)
+ sigma_spec_lm(nd,lth,kr) = sigma_spec_lm(nd,lth,kr) &
+ & + half * (tmp_lm + spectr_pre_l(nd,lth,kr)) &
+ & * (time_sph + pre_time)
+!
+ spectr_pre_l(nd,lth,kr) = tmp_l
+ spectr_pre_m(nd,lth,kr) = tmp_m
+ spectr_pre_lm(nd,lth,kr) = tmp_lm
+ end do
+!$omp end do nowait
+ end do
+ end do
+!$omp end parallel
+!
+ pre_time = time_sph
+!
+ end subroutine sum_deviation_ene_sph
+!
+! --------------------------------------------------------------------
+!
+ subroutine divide_average_ene_sph
+!
+ integer(kind = kint) :: kr, nd, lth
+!
+!
+!$omp parallel private(kr,lth,nd)
+ do kr = 1, nri_sph
+!$omp do
+ do nd = 1, ncomp_sph_spec
+ ave_spec_t(nd,kr) = ave_spec_t(nd,kr) / (time_sph - time_ini)
+ end do
+!$omp end do nowait
+!
+ do lth = 0, ltr_sph
+!$omp do
+ do nd = 1, ncomp_sph_spec
+ ave_spec_l(nd,lth,kr) = ave_spec_l(nd,lth,kr) &
+ & / (time_sph - time_ini)
+ ave_spec_m(nd,lth,kr) = ave_spec_m(nd,lth,kr) &
+ & / (time_sph - time_ini)
+ ave_spec_lm(nd,lth,kr) = ave_spec_lm(nd,lth,kr) &
+ & / (time_sph - time_ini)
+ end do
+!$omp end do nowait
+ end do
+ end do
+!$omp end parallel
+!
+ end subroutine divide_average_ene_sph
+!
+! --------------------------------------------------------------------
+!
+ subroutine divide_deviation_ene_sph
+!
+ integer(kind = kint) :: kr, nd, lth
+!
+!
+!$omp parallel private(kr,lth,nd)
+ do kr = 1, nri_sph
+!$omp do
+ do nd = 1, ncomp_sph_spec
+ sigma_spec_t(nd,kr) = sqrt(sigma_spec_t(nd,kr) &
+ & / (time_sph - time_ini))
+ end do
+!$omp end do nowait
+!
+ do lth = 0, ltr_sph
+!$omp do
+ do nd = 1, ncomp_sph_spec
+ sigma_spec_l(nd,lth,kr) = sqrt(sigma_spec_l(nd,lth,kr) &
+ & / (time_sph - time_ini))
+ sigma_spec_m(nd,lth,kr) = sqrt(sigma_spec_m(nd,lth,kr) &
+ & / (time_sph - time_ini))
+ sigma_spec_lm(nd,lth,kr) = sqrt(sigma_spec_lm(nd,lth,kr) &
+ & / (time_sph - time_ini))
+ end do
+!$omp end do nowait
+ end do
+ end do
+!$omp end parallel
+!
+ end subroutine divide_deviation_ene_sph
+!
+! --------------------------------------------------------------------
+!
+ subroutine output_tave_ene_sph_data
+!
+!
+ call open_tave_ene_spec_data
+!
+ if(iflag_sph_ene_file .eq. 1) then
+ call write_tave_vol_sph_data(nri_sph, ltr_sph, ncomp_sph_spec, &
+ & ave_spec_t, ave_spec_l, ave_spec_m, ave_spec_lm)
+ else
+ call write_tave_layer_sph_data(nri_sph, ltr_sph, &
+ & ncomp_sph_spec, ave_spec_t, ave_spec_l, ave_spec_m, &
+ & ave_spec_lm)
+ end if
+!
+ call close_ene_spec_data
+!
+ end subroutine output_tave_ene_sph_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine output_tsigma_ene_sph_data
+!
+!
+ call open_tsigma_ene_spec_data
+!
+ if(iflag_sph_ene_file .eq. 1) then
+ call write_tave_vol_sph_data(nri_sph, ltr_sph, ncomp_sph_spec, &
+ & sigma_spec_t, sigma_spec_l, sigma_spec_m, sigma_spec_lm)
+ else
+ call write_tave_layer_sph_data(nri_sph, ltr_sph, &
+ & ncomp_sph_spec, sigma_spec_t, sigma_spec_l, sigma_spec_m, &
+ & sigma_spec_lm)
+ end if
+!
+ call close_ene_spec_data
+!
+ end subroutine output_tsigma_ene_sph_data
+!
+! --------------------------------------------------------------------
+!
+ end module cal_tave_sph_ene_spectr
diff --git a/src/programs/data_utilities/TIME_HISTORIES/m_sph_ene_spectra.f90 b/src/programs/data_utilities/TIME_HISTORIES/m_sph_ene_spectra.f90
new file mode 100644
index 0000000..c84354a
--- /dev/null
+++ b/src/programs/data_utilities/TIME_HISTORIES/m_sph_ene_spectra.f90
@@ -0,0 +1,442 @@
+!>@file m_sph_ene_spectra.f90
+!! module m_sph_ene_spectra
+!!
+!! @author H. Matsui
+!! @date Programmed in Nov., 2007
+!!
+!
+!> @brief Time average spherical harmonics spectrum data
+!!
+!!@verbatim
+!! subroutine allocate_sph_espec_name
+!! subroutine allocate_tave_sph_espec_data
+!! subroutine deallocate_tave_sph_espec_data
+!!
+!! subroutine select_sph_ene_spec_data_file
+!! subroutine set_org_ene_spec_file_name
+!!
+!! subroutine open_org_ene_spec_data
+!! subroutine close_ene_spec_data
+!!
+!! subroutine read_org_layer_ene_data(istep, ierr)
+!! subroutine read_org_volume_ene_data(istep, ierr)
+!!@endverbatim
+!
+ module m_sph_ene_spectra
+!
+ use m_precision
+ use m_constants
+!
+ implicit none
+!
+!
+ integer(kind = kint) :: ltr_sph, nri_sph
+ integer(kind = kint) :: ncomp_sph_spec, num_time_labels
+ character(len = kchara), allocatable :: ene_sph_spec_name(:)
+!
+ integer(kind = kint) :: istep_st, istep_ed, istep_read
+ integer(kind = kint) :: ist_true, ied_true
+ real(kind = kreal) :: time_sph, time_ini
+!
+ integer(kind = kint), allocatable :: kr_sph(:)
+ real(kind = kreal), allocatable :: spectr_t(:,:)
+ real(kind = kreal), allocatable :: spectr_l(:,:,:)
+ real(kind = kreal), allocatable :: spectr_m(:,:,:)
+ real(kind = kreal), allocatable :: spectr_lm(:,:,:)
+!
+ real(kind = kreal), allocatable :: spectr_pre_t(:,:)
+ real(kind = kreal), allocatable :: spectr_pre_l(:,:,:)
+ real(kind = kreal), allocatable :: spectr_pre_m(:,:,:)
+ real(kind = kreal), allocatable :: spectr_pre_lm(:,:,:)
+!
+ real(kind = kreal) :: pre_time
+!
+ integer(kind = kint) :: iflag_sph_ene_file
+ integer(kind = kint) :: ilayer_sph_ene
+!
+! data file ID
+!
+ integer(kind = kint), parameter :: id_file_rms_l = 31
+ integer(kind = kint), parameter :: id_file_rms_m = 32
+ integer(kind = kint), parameter :: id_file_rms_lm = 33
+ integer(kind = kint), parameter :: id_file_rms = 34
+!
+! data file name
+!
+ character(len = kchara) :: fhead_rms_vol = 'sph_pwr_volume'
+ character(len = kchara) :: fhead_rms_layer = 'sph_pwr_layer'
+!
+ character(len = kchara) :: fname_org_rms_l
+ character(len = kchara) :: fname_org_rms_m
+ character(len = kchara) :: fname_org_rms_lm
+ character(len = kchara) :: fname_org_rms
+!
+ private :: fhead_rms_vol, fhead_rms_layer
+ private :: ilayer_sph_ene
+ private :: read_ene_spectr_header
+!
+! --------------------------------------------------------------------
+!
+ contains
+!
+! --------------------------------------------------------------------
+!
+ subroutine allocate_sph_espec_name
+!
+!
+ allocate( ene_sph_spec_name(ncomp_sph_spec+num_time_labels) )
+!
+ end subroutine allocate_sph_espec_name
+!
+! --------------------------------------------------------------------
+!
+ subroutine allocate_sph_espec_data
+!
+!
+ allocate( kr_sph(nri_sph) )
+!
+ allocate( spectr_t(ncomp_sph_spec,nri_sph) )
+ allocate( spectr_l(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+ allocate( spectr_m(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+ allocate( spectr_lm(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+!
+ allocate( spectr_pre_t(ncomp_sph_spec,nri_sph) )
+ allocate( spectr_pre_l(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+ allocate( spectr_pre_m(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+ allocate( spectr_pre_lm(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+!
+ spectr_t = zero
+ spectr_l = zero
+ spectr_m = zero
+ spectr_lm = zero
+!
+ spectr_pre_t = zero
+ spectr_pre_l = zero
+ spectr_pre_m = zero
+ spectr_pre_lm = zero
+!
+ end subroutine allocate_sph_espec_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine deallocate_tave_sph_espec_data
+!
+ deallocate(ene_sph_spec_name, kr_sph)
+ deallocate(spectr_t, spectr_l, spectr_m, spectr_lm)
+ deallocate(spectr_pre_t, spectr_pre_l)
+ deallocate(spectr_pre_m, spectr_pre_lm)
+!
+ end subroutine deallocate_tave_sph_espec_data
+!
+! --------------------------------------------------------------------
+! --------------------------------------------------------------------
+!
+ subroutine select_sph_ene_spec_data_file
+!
+!
+ write(*,*) ' Choose data type to taking average'
+ write(*,*) ' 1: volume average spectrum '
+ write(*,*) ' 2: spectrum on each layer '
+ write(*,*) ' 3: spectrum on specific layer '
+ read(*,*) iflag_sph_ene_file
+!
+ write(*,*) 'enter file header for averaging'
+ if (iflag_sph_ene_file .eq. 1) then
+ read(*,*) fhead_rms_vol
+ else
+ read(*,*) fhead_rms_layer
+ end if
+!
+ if (iflag_sph_ene_file .eq. 3) then
+ write(*,*) ' Choose layer number'
+ read(*,*) ilayer_sph_ene
+ end if
+!
+ end subroutine select_sph_ene_spec_data_file
+!
+! --------------------------------------------------------------------
+!
+ subroutine set_org_ene_spec_file_name
+!
+ use set_parallel_file_name
+!
+ character(len = kchara) :: fname_tmp
+!
+ if (iflag_sph_ene_file .eq. 1) then
+ write(fname_org_rms_l, '(a,a6)') &
+ & trim(fhead_rms_vol), '_l.dat'
+ write(fname_org_rms_m, '(a,a6)') &
+ & trim(fhead_rms_vol), '_m.dat'
+ write(fname_org_rms_lm,'(a,a7)') &
+ & trim(fhead_rms_vol), '_lm.dat'
+ call add_dat_extension(fhead_rms_vol, fname_org_rms)
+ else if (iflag_sph_ene_file .eq. 2) then
+ write(fname_org_rms_l, '(a,a6)') &
+ & trim(fhead_rms_layer), '_l.dat'
+ write(fname_org_rms_m, '(a,a6)') &
+ & trim(fhead_rms_layer), '_m.dat'
+ write(fname_org_rms_lm,'(a,a7)') &
+ & trim(fhead_rms_layer), '_lm.dat'
+ call add_dat_extension(fhead_rms_layer, fname_org_rms)
+ else if (iflag_sph_ene_file .eq. 3) then
+ write(fname_org_rms_l, '(a,a2)') fhead_rms_layer, '_l'
+ call add_int_suffix(ilayer_sph_ene, &
+ & fname_org_rms_l, fname_tmp)
+ call add_dat_extension(fname_tmp, fname_org_rms_l)
+!
+ write(fname_org_rms_m, '(a,a2)') fhead_rms_layer, '_m'
+ call add_int_suffix(ilayer_sph_ene, &
+ & fname_org_rms_m, fname_tmp)
+ call add_dat_extension(fname_tmp, fname_org_rms_m)
+!
+ write(fname_org_rms_lm, '(a,a3)') fhead_rms_layer, '_lm'
+ call add_int_suffix(ilayer_sph_ene, &
+ & fname_org_rms_lm, fname_tmp)
+ call add_dat_extension(fname_tmp, fname_org_rms_lm)
+!
+ write(fname_org_rms, '(a,a3)') fhead_rms_layer, '_lm'
+ call add_int_suffix(ilayer_sph_ene, &
+ & fname_org_rms, fname_tmp)
+ call add_dat_extension(fname_tmp, fname_org_rms)
+ end if
+!
+ end subroutine set_org_ene_spec_file_name
+!
+! --------------------------------------------------------------------
+!
+ subroutine open_org_ene_spec_data
+!
+!
+ open(id_file_rms, file=fname_org_rms)
+ open(id_file_rms_l, file=fname_org_rms_l)
+ open(id_file_rms_m, file=fname_org_rms_m)
+ open(id_file_rms_lm,file=fname_org_rms_lm)
+!
+ call read_ene_spectr_header(id_file_rms, ione)
+ call read_ene_spectr_header(id_file_rms_l, izero)
+ call read_ene_spectr_header(id_file_rms_m, izero)
+ call read_ene_spectr_header(id_file_rms_lm, izero)
+!
+ end subroutine open_org_ene_spec_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine close_ene_spec_data
+!
+!
+ close(id_file_rms)
+ close(id_file_rms_l)
+ close(id_file_rms_m)
+ close(id_file_rms_lm)
+!
+ end subroutine close_ene_spec_data
+!
+! --------------------------------------------------------------------
+! --------------------------------------------------------------------
+!
+ subroutine read_ene_spectr_header(id_file, iflag_total)
+!
+ use skip_comment_f
+!
+ integer(kind = kint), intent(in) :: id_file, iflag_total
+ character(len=255) :: character_4_read
+ integer(kind = kint) :: num, itmp, nfld, i
+!
+!
+ call skip_comment(character_4_read, id_file)
+ call skip_comment(character_4_read, id_file)
+ call skip_comment(character_4_read, id_file)
+ call skip_comment(character_4_read, id_file)
+ call skip_comment(character_4_read, id_file)
+ read(id_file,*) nfld, ncomp_sph_spec
+ read(id_file,*) (itmp,i=1,nfld)
+!
+ num = size(ene_sph_spec_name) - iflag_total
+ read(id_file,*) ene_sph_spec_name(1:num)
+!
+ end subroutine read_ene_spectr_header
+!
+! --------------------------------------------------------------------
+! --------------------------------------------------------------------
+!
+ subroutine read_org_layer_ene_data(istep, ierr)
+!
+ integer(kind = kint), intent(inout) :: istep, ierr
+ integer(kind = kint) :: itmp
+ integer(kind = kint) :: kr, lth
+!
+!
+ ierr = 0
+ do kr = 1, nri_sph
+ read(id_file_rms,*,err=99,end=99) istep, time_sph, &
+ & kr_sph(kr), spectr_t(1:ncomp_sph_spec,kr)
+ do lth = 0, ltr_sph
+ read(id_file_rms_l,*,err=99,end=99) istep, time_sph, &
+ & itmp, itmp, spectr_l(1:ncomp_sph_spec,lth,kr)
+ read(id_file_rms_m,*,err=99,end=99) istep, time_sph, &
+ & itmp, itmp, spectr_m(1:ncomp_sph_spec,lth,kr)
+ read(id_file_rms_lm,*,err=99,end=99) istep, time_sph, &
+ & itmp, itmp, spectr_lm(1:ncomp_sph_spec,lth,kr)
+ end do
+ end do
+ return
+!
+ 99 continue
+ ierr = 1
+ return
+!
+ end subroutine read_org_layer_ene_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine read_org_volume_ene_data(istep, ierr)
+!
+ integer(kind = kint), intent(inout) :: istep, ierr
+ integer(kind = kint) :: itmp
+ integer(kind = kint) :: lth
+!
+!
+ ierr = 0
+ read(id_file_rms,*,err=99,end=99) istep, time_sph, &
+ & spectr_t(1:ncomp_sph_spec,ione)
+ do lth = 0, ltr_sph
+ read(id_file_rms_l,*,err=99,end=99) istep, time_sph, &
+ & itmp, spectr_l(1:ncomp_sph_spec,lth,ione)
+ read(id_file_rms_m,*,err=99,end=99) istep, time_sph, &
+ & itmp, spectr_m(1:ncomp_sph_spec,lth,ione)
+ read(id_file_rms_lm,*,err=99,end=99) istep, time_sph, &
+ & itmp, spectr_lm(1:ncomp_sph_spec,lth,ione)
+ end do
+ return
+!
+ 99 continue
+ ierr = 1
+ return
+!
+ end subroutine read_org_volume_ene_data
+!
+! --------------------------------------------------------------------
+! --------------------------------------------------------------------
+!
+ subroutine count_degree_on_layer_data
+!
+ use skip_comment_f
+!
+ character(len=255) :: character_4_read
+ integer(kind = kint) :: num, itmp, nfld, i
+!
+!
+!
+ open (id_file_rms_l,file=fname_org_rms_l)
+!
+ call skip_comment(character_4_read, id_file_rms_l)
+ read(id_file_rms_l,*) nri_sph, ltr_sph
+ write(*,*) 'ltr_sph', ltr_sph
+ write(*,*) 'nri_sph', nri_sph
+ call skip_comment(character_4_read, id_file_rms_l)
+ call skip_comment(character_4_read, id_file_rms_l)
+ call skip_comment(character_4_read, id_file_rms_l)
+ read(id_file_rms_l,*) nfld, ncomp_sph_spec
+!
+ num_time_labels = 4
+ write(*,*) 'ncomp_sph_spec', ncomp_sph_spec
+ call allocate_sph_espec_name
+!
+!
+ read(id_file_rms_l,*) (itmp,i=1,nfld)
+!
+ num = size(ene_sph_spec_name)
+ write(*,*) 'num vector', num
+ read(id_file_rms_l,*) ene_sph_spec_name(1:num)
+ do i = 1, NUM
+ write(*,*) i, trim(ene_sph_spec_name(i))
+ end do
+!
+ close(id_file_rms_l)
+!
+ end subroutine count_degree_on_layer_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine count_degree_one_layer_data
+!
+ use skip_comment_f
+!
+ character(len=255) :: character_4_read
+ integer(kind = kint) :: num, itmp, nfld, i
+!
+!
+!
+ open (id_file_rms_l,file=fname_org_rms_l)
+!
+ call skip_comment(character_4_read, id_file_rms_l)
+ read(id_file_rms_l,*) nri_sph, ltr_sph
+ nri_sph = 1
+ write(*,*) 'ltr_sph', ltr_sph
+ write(*,*) 'nri_sph', nri_sph
+ call skip_comment(character_4_read, id_file_rms_l)
+ call skip_comment(character_4_read, id_file_rms_l)
+ call skip_comment(character_4_read, id_file_rms_l)
+ read(id_file_rms_l,*) nfld, ncomp_sph_spec
+!
+ num_time_labels = 4
+ write(*,*) 'ncomp_sph_spec', ncomp_sph_spec
+ call allocate_sph_espec_name
+!
+!
+ read(id_file_rms_l,*) (itmp,i=1,nfld)
+!
+ num = size(ene_sph_spec_name)
+ write(*,*) 'num vector', num
+ read(id_file_rms_l,*) ene_sph_spec_name(1:num)
+ do i = 1, NUM
+ write(*,*) i, trim(ene_sph_spec_name(i))
+ end do
+!
+ close(id_file_rms_l)
+!
+ end subroutine count_degree_one_layer_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine count_degree_on_volume_data
+!
+ use skip_comment_f
+!
+ character(len=255) :: character_4_read
+ integer(kind = kint) :: num, itmp, nfld, i
+!
+!
+ open (id_file_rms_l,file=fname_org_rms_l)
+!
+ call skip_comment(character_4_read, id_file_rms_l)
+ read(id_file_rms_l,*) itmp, ltr_sph
+ write(*,*) 'ltr_sph', ltr_sph
+ call skip_comment(character_4_read, id_file_rms_l)
+ call skip_comment(character_4_read, id_file_rms_l)
+ call skip_comment(character_4_read, id_file_rms_l)
+ read(id_file_rms_l,*) nfld, ncomp_sph_spec
+!
+ num_time_labels = 3
+ nri_sph = 1
+ write(*,*) 'ncomp_sph_spec', ncomp_sph_spec
+ write(*,*) 'nri_sph', nri_sph
+ call allocate_sph_espec_name
+!
+!
+ read(id_file_rms_l,*) (itmp,i=1,nfld)
+!
+ num = size(ene_sph_spec_name)
+ write(*,*) 'num v', num
+ read(id_file_rms_l,*) ene_sph_spec_name(1:num)
+ DO I = 1, NUM
+ write(*,*) i, trim(ene_sph_spec_name(i))
+ end do
+!
+ close(id_file_rms_l)
+!
+ end subroutine count_degree_on_volume_data
+!
+! --------------------------------------------------------------------
+!
+ end module m_sph_ene_spectra
diff --git a/src/programs/data_utilities/TIME_HISTORIES/m_tave_sph_ene_spectr.f90 b/src/programs/data_utilities/TIME_HISTORIES/m_tave_sph_ene_spectr.f90
new file mode 100644
index 0000000..ea13a17
--- /dev/null
+++ b/src/programs/data_utilities/TIME_HISTORIES/m_tave_sph_ene_spectr.f90
@@ -0,0 +1,280 @@
+!>@file m_tave_sph_ene_spectr.f90
+!! module m_tave_sph_ene_spectr
+!!
+!! @author H. Matsui
+!! @date Programmed in Nov., 2007
+!!
+!
+!> @brief Time average spherical harmonics spectrum data
+!!
+!!@verbatim
+!! subroutine allocate_tave_sph_espec_data
+!! subroutine deallocate_tave_sph_espec_data
+!!
+!! subroutine reset_tave_sph_espec_data
+!! subroutine reset_tsigma_sph_espec_data
+!!
+!! subroutine open_tave_ene_spec_data
+!! subroutine open_tsigma_ene_spec_data
+!!
+!! subroutine write_tave_vol_sph_data(nri, ltr, ncomp, &
+!! & spec, spec_l, spec_m, spec_lm)
+!! subroutine write_tave_layer_sph_data(nri, ltr, ncomp, &
+!! & spec, spec_l, spec_m, spec_lm)
+!!@endverbatim
+!
+ module m_tave_sph_ene_spectr
+!
+ use m_precision
+ use m_constants
+!
+ implicit none
+!
+!
+ real(kind = kreal), allocatable :: ave_spec_t(:,:)
+ real(kind = kreal), allocatable :: ave_spec_l(:,:,:)
+ real(kind = kreal), allocatable :: ave_spec_m(:,:,:)
+ real(kind = kreal), allocatable :: ave_spec_lm(:,:,:)
+!
+ real(kind = kreal), allocatable :: sigma_spec_t(:,:)
+ real(kind = kreal), allocatable :: sigma_spec_l(:,:,:)
+ real(kind = kreal), allocatable :: sigma_spec_m(:,:,:)
+ real(kind = kreal), allocatable :: sigma_spec_lm(:,:,:)
+!
+! data file ID
+!
+ integer(kind = kint), parameter :: id_tave_rms_l = 31
+ integer(kind = kint), parameter :: id_tave_rms_m = 32
+ integer(kind = kint), parameter :: id_tave_rms_lm = 33
+ integer(kind = kint), parameter :: id_tave_rms = 34
+!
+ private :: write_average_ene_sph_head
+!
+! --------------------------------------------------------------------
+!
+ contains
+!
+! --------------------------------------------------------------------
+!
+ subroutine allocate_tave_sph_espec_data
+!
+ use m_sph_ene_spectra
+!
+ allocate( ave_spec_t(ncomp_sph_spec,nri_sph) )
+ allocate( ave_spec_l(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+ allocate( ave_spec_m(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+ allocate( ave_spec_lm(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+!
+ allocate( sigma_spec_t(ncomp_sph_spec,nri_sph) )
+ allocate( sigma_spec_l(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+ allocate( sigma_spec_m(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+ allocate( sigma_spec_lm(ncomp_sph_spec,0:ltr_sph,nri_sph) )
+!
+ ave_spec_t = 0.0d0
+ ave_spec_l = 0.0d0
+ ave_spec_m = 0.0d0
+ ave_spec_lm = 0.0d0
+!
+ sigma_spec_t = 0.0d0
+ sigma_spec_l = 0.0d0
+ sigma_spec_m = 0.0d0
+ sigma_spec_lm = 0.0d0
+!
+ end subroutine allocate_tave_sph_espec_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine deallocate_tave_sph_espec_data
+!
+ deallocate(ave_spec_t, ave_spec_l, ave_spec_m, ave_spec_lm)
+ deallocate(sigma_spec_l, sigma_spec_m)
+ deallocate(sigma_spec_t, sigma_spec_lm)
+!
+ end subroutine deallocate_tave_sph_espec_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine reset_tave_sph_espec_data
+!
+ use m_sph_ene_spectra
+!
+!
+ sigma_spec_t = 0.0d0
+ sigma_spec_l = 0.0d0
+ sigma_spec_m = 0.0d0
+ sigma_spec_lm = 0.0d0
+!
+ end subroutine reset_tave_sph_espec_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine reset_tsigma_sph_espec_data
+!
+!
+ sigma_spec_t = 0.0d0
+ sigma_spec_l = 0.0d0
+ sigma_spec_m = 0.0d0
+ sigma_spec_lm = 0.0d0
+!
+ end subroutine reset_tsigma_sph_espec_data
+!
+! --------------------------------------------------------------------
+! --------------------------------------------------------------------
+!
+ subroutine open_tave_ene_spec_data
+!
+ use m_sph_ene_spectra
+!
+ character(len = kchara) :: fname_tave_rms_l
+ character(len = kchara) :: fname_tave_rms_m
+ character(len = kchara) :: fname_tave_rms_lm
+ character(len = kchara) :: fname_tave_rms
+!
+!
+ write(fname_tave_rms_l, '(a6,a)') &
+ & 't_ave_', trim(fname_org_rms_l)
+ write(fname_tave_rms_m, '(a6,a)') &
+ & 't_ave_', trim(fname_org_rms_m)
+ write(fname_tave_rms_lm,'(a6,a)') &
+ & 't_ave_', trim(fname_org_rms_lm)
+ write(fname_tave_rms, '(a6,a)') &
+ & 't_ave_', trim(fname_org_rms)
+!
+ open(id_tave_rms, file=fname_tave_rms)
+ open(id_tave_rms_l, file=fname_tave_rms_l)
+ open(id_tave_rms_m, file=fname_tave_rms_m)
+ open(id_tave_rms_lm,file=fname_tave_rms_lm)
+!
+ call write_average_ene_sph_head
+!
+ end subroutine open_tave_ene_spec_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine open_tsigma_ene_spec_data
+!
+ use m_sph_ene_spectra
+!
+ character(len = kchara) :: fname_tave_rms_l
+ character(len = kchara) :: fname_tave_rms_m
+ character(len = kchara) :: fname_tave_rms_lm
+ character(len = kchara) :: fname_tave_rms
+!
+!
+ write(fname_tave_rms_l, '(a8,a)') &
+ & 't_sigma_', trim(fname_org_rms_l)
+ write(fname_tave_rms_m, '(a8,a)') &
+ & 't_sigma_', trim(fname_org_rms_m)
+ write(fname_tave_rms_lm,'(a8,a)') &
+ & 't_sigma_', trim(fname_org_rms_lm)
+ write(fname_tave_rms, '(a8,a)') &
+ & 't_sigma_', trim(fname_org_rms)
+!
+ open(id_tave_rms, file=fname_tave_rms)
+ open(id_tave_rms_l, file=fname_tave_rms_l)
+ open(id_tave_rms_m, file=fname_tave_rms_m)
+ open(id_tave_rms_lm,file=fname_tave_rms_lm)
+!
+ call write_average_ene_sph_head
+!
+ end subroutine open_tsigma_ene_spec_data
+!
+! --------------------------------------------------------------------
+! --------------------------------------------------------------------
+!
+ subroutine write_average_ene_sph_head
+!
+ use m_sph_ene_spectra
+ use write_field_labels
+!
+ integer(kind = kint) :: num
+!
+!
+ num = ncomp_sph_spec + num_time_labels
+ write(ene_sph_spec_name(num_time_labels),'(a)') 'degree'
+ call write_multi_labels(id_tave_rms, num, ene_sph_spec_name)
+ call write_multi_labels(id_tave_rms_l, num, ene_sph_spec_name)
+!
+ write(ene_sph_spec_name(num_time_labels),'(a)') 'order'
+ call write_multi_labels(id_tave_rms_m, num, ene_sph_spec_name)
+!
+ write(ene_sph_spec_name(num_time_labels),'(a)') &
+ & 'diff_deg_order'
+ call write_multi_labels(id_tave_rms_lm, num, ene_sph_spec_name)
+!
+ write(id_tave_rms,* )
+ write(id_tave_rms_l,* )
+ write(id_tave_rms_m,* )
+ write(id_tave_rms_lm,*)
+!
+ end subroutine write_average_ene_sph_head
+!
+! --------------------------------------------------------------------
+!
+ subroutine write_tave_vol_sph_data(nri, ltr, ncomp, &
+ & spec, spec_l, spec_m, spec_lm)
+!
+ use m_sph_ene_spectra
+!
+ integer(kind = kint), intent(in) :: nri, ltr, ncomp
+ real(kind = kreal), intent(inout) :: spec(ncomp,nri)
+ real(kind = kreal), intent(inout) :: spec_l(ncomp,0:ltr,nri)
+ real(kind = kreal), intent(inout) :: spec_m(ncomp,0:ltr,nri)
+ real(kind = kreal), intent(inout) :: spec_lm(ncomp,0:ltr,nri)
+!
+ integer(kind = kint) :: lth
+!
+!
+ write(id_tave_rms,'(i10,1pE25.15e3,i10,1p255E25.15e3)') &
+ & ied_true, time_sph, izero, spec(1:ncomp,1)
+!
+ do lth = 0, ltr
+ write(id_tave_rms_l,'(i10,1pE25.15e3,i10,1p255E25.15e3)') &
+ & ied_true, time_sph, lth, spec_l(1:ncomp,lth,1)
+ write(id_tave_rms_m,'(i10,1pE25.15e3,i10,1p255E25.15e3)') &
+ & ied_true, time_sph, lth, spec_m(1:ncomp,lth,1)
+ write(id_tave_rms_lm,'(i10,1pE25.15e3,i10,1p255E25.15e3)') &
+ & ied_true, time_sph, lth, spec_lm(1:ncomp,lth,1)
+ end do
+!
+ end subroutine write_tave_vol_sph_data
+!
+! --------------------------------------------------------------------
+!
+ subroutine write_tave_layer_sph_data(nri, ltr, ncomp, &
+ & spec, spec_l, spec_m, spec_lm)
+!
+ use m_sph_ene_spectra
+!
+ integer(kind = kint), intent(in) :: nri, ltr, ncomp
+ real(kind = kreal), intent(inout) :: spec(ncomp,nri)
+ real(kind = kreal), intent(inout) :: spec_l(ncomp,0:ltr,nri)
+ real(kind = kreal), intent(inout) :: spec_m(ncomp,0:ltr,nri)
+ real(kind = kreal), intent(inout) :: spec_lm(ncomp,0:ltr,nri)
+!
+ integer(kind = kint) :: kr, lth
+!
+!
+ do kr = 1, nri
+ write(id_tave_rms,'(i10,1pE25.15e3,2i10,1p255E25.15e3)') &
+ & ied_true, time_sph, kr_sph(kr), izero, &
+ & spec(1:ncomp,kr)
+!
+ do lth = 0, ltr
+ write(id_tave_rms_l,'(i10,1pE25.15e3,2i10,1p255E25.15e3)') &
+ & ied_true, time_sph, kr_sph(kr), lth, &
+ & spec_l(1:ncomp,lth,kr)
+ write(id_tave_rms_m,'(i10,1pE25.15e3,2i10,1p255E25.15e3)') &
+ & ied_true, time_sph, kr_sph(kr), lth, &
+ & spec_m(1:ncomp,lth,kr)
+ write(id_tave_rms_lm,'(i10,1pE25.15e3,2i10,1p255E25.15e3)') &
+ & ied_true, time_sph, kr_sph(kr), lth, &
+ & spec_lm(1:ncomp,lth,kr)
+ end do
+ end do
+!
+ end subroutine write_tave_layer_sph_data
+!
+! --------------------------------------------------------------------
+!
+ end module m_tave_sph_ene_spectr
diff --git a/src/programs/data_utilities/TIME_HISTORIES/t_average_sph_ene_spec.f90 b/src/programs/data_utilities/TIME_HISTORIES/t_average_sph_ene_spec.f90
new file mode 100644
index 0000000..bc7694b
--- /dev/null
+++ b/src/programs/data_utilities/TIME_HISTORIES/t_average_sph_ene_spec.f90
@@ -0,0 +1,123 @@
+!>@file t_average_sph_ene_spec.f90
+!! program t_average_sph_ene_spec
+!!
+!! @author H. Matsui
+!! @date Programmed in Nov., 2007
+!!
+!
+!> @brief Evaluate time average and standard deviation
+!! from spherical harmonics spectrum data
+!
+ program t_average_sph_ene_spec
+!
+ use m_precision
+!
+ use m_tave_sph_ene_spectr
+ use cal_tave_sph_ene_spectr
+!
+ implicit none
+!
+!
+ integer(kind = kint) :: ierr
+ integer(kind = kint) :: icou, istep
+ real(kind = kreal) :: start_time, end_time
+!
+!
+ call select_sph_ene_spec_data_file
+ call set_org_ene_spec_file_name
+!
+ write(*,*) 'Input start and end time'
+ read(*,*) start_time, end_time
+!
+ if(iflag_sph_ene_file .eq. 1) then
+ call count_degree_on_volume_data
+ else if(iflag_sph_ene_file .eq. 2) then
+ call count_degree_on_layer_data
+ else
+ call count_degree_one_layer_data
+ end if
+ call allocate_sph_espec_data
+ call allocate_tave_sph_espec_data
+!
+! Evaluate time average
+!
+ call open_org_ene_spec_data
+!
+ ist_true = -1
+ icou = 0
+ do
+ if(iflag_sph_ene_file .eq. 1) then
+ call read_org_volume_ene_data(istep, ierr)
+ else
+ call read_org_layer_ene_data(istep, ierr)
+ end if
+ if(ierr.gt.0) go to 99
+!
+ if (time_sph .ge. start_time) then
+ if (ist_true .eq. -1) then
+ ist_true = istep
+ end if
+ icou = icou + 1
+ ied_true = istep
+!
+ call sum_average_ene_sph
+ if (icou .eq. 1) then
+ time_ini = time_sph
+ call reset_tave_sph_espec_data
+ end if
+ end if
+!
+ if (time_sph .ge. end_time) exit
+!
+ write(*,*) 'step', istep, 'averaging finished. Count: ', icou
+ end do
+ 99 continue
+ call close_ene_spec_data
+!
+ call divide_average_ene_sph
+ call output_tave_ene_sph_data
+!
+ call close_ene_spec_data
+!
+! Evaluate standard deviation
+!
+ call open_org_ene_spec_data
+!
+ ist_true = -1
+ icou = 0
+ do
+ if(iflag_sph_ene_file .eq. 1) then
+ call read_org_volume_ene_data(istep, ierr)
+ else
+ call read_org_layer_ene_data(istep, ierr)
+ end if
+ if(ierr.gt.0) go to 98
+!
+ if (time_sph .ge. start_time) then
+ if (ist_true .eq. -1) then
+ ist_true = istep
+ end if
+ icou = icou + 1
+ ied_true = istep
+!
+ call sum_deviation_ene_sph
+ if (icou .eq. 1) then
+ time_ini = time_sph
+ call reset_tsigma_sph_espec_data
+ end if
+ end if
+!
+ if (time_sph .ge. end_time) exit
+!
+ write(*,*) 'step', istep, 'deviation finished. Count: ', icou
+!
+ end do
+ 98 continue
+ call close_ene_spec_data
+!
+ call divide_deviation_ene_sph
+ call output_tsigma_ene_sph_data
+
+!
+ stop
+ end program t_average_sph_ene_spec
diff --git a/src/programs/data_utilities/TIME_HISTORIES/tave_picked_sph_spec_data.f90 b/src/programs/data_utilities/TIME_HISTORIES/tave_picked_sph_spec_data.f90
new file mode 100644
index 0000000..f28d0d9
--- /dev/null
+++ b/src/programs/data_utilities/TIME_HISTORIES/tave_picked_sph_spec_data.f90
@@ -0,0 +1,171 @@
+!tave_picked_sph_spec_data.f90
+! program tave_picked_sph_spec_data
+!
+! programmed by H.Matsui on Dec., 2012
+!
+ program tave_picked_sph_spec_data
+!
+ use m_precision
+ use m_constants
+!
+ use m_pickup_sph_spectr_data
+!
+ implicit none
+!
+ real(kind = kreal), allocatable :: prev_spec(:,:)
+ real(kind = kreal), allocatable :: ave_spec(:,:)
+ real(kind = kreal), allocatable :: sdev_spec(:,:)
+ character(len=kchara) :: evo_header
+ character(len=kchara) :: tave_header
+ character(len=kchara) :: sdev_header
+ integer(kind = kint), parameter :: id_pick = 15
+!
+ integer(kind = kint) :: i_step, ierr, num, icou, ipick, nd
+ integer(kind = kint) :: istep_start
+ real(kind = kreal) :: acou, time, prev_time
+ real(kind = kint) :: start_time, end_time, true_start
+!
+!
+ write(*,*) 'Input picked spectr evolution file header'
+ read(5,*) evo_header
+!
+ pickup_sph_head = evo_header
+ write(tave_header,'(a6,a)') 't_ave_', trim(evo_header)
+ write(sdev_header,'(a8,a)') 't_sigma_', trim(evo_header)
+!
+ write(*,*) 'Input start and end time'
+ read(5,*) start_time, end_time
+!
+ call open_sph_spec_read_monitor(id_pick)
+!
+ num = ntot_pick_sph_mode*num_pick_layer
+ allocate( prev_spec(ntot_comp_pick_sph,num) )
+ allocate( ave_spec(ntot_comp_pick_sph,num) )
+ allocate( sdev_spec(ntot_comp_pick_sph,num) )
+ prev_spec = 0.0d0
+ ave_spec = 0.0d0
+ sdev_spec = 0.0d0
+!
+! Evaluate time average
+!
+ icou = 0
+ do
+ call read_sph_spec_4_monitor(id_pick, i_step, time, ierr)
+ if(ierr .gt. 0) exit
+!
+ if(time .ge. start_time) then
+ do ipick = 1, num_pick_sph_mode*num_pick_layer
+ do nd = 1, ntot_comp_pick_sph
+ prev_spec(nd,ipick) = d_rj_pick_sph_gl(nd,ipick)
+ end do
+ end do
+!
+ if(icou .eq. 0) then
+ true_start = time
+ else
+ do ipick = 1, num_pick_sph_mode*num_pick_layer
+ do nd = 1, ntot_comp_pick_sph
+ ave_spec(nd,ipick) = ave_spec(nd,ipick) + half &
+ & * (d_rj_pick_sph_gl(nd,ipick) + prev_spec(nd,ipick)) &
+ & * (time - prev_time)
+ end do
+ end do
+ end if
+!
+ icou = icou + 1
+ write(*,*) 'step ', i_step, &
+ & ' is added for time average: count is ', icou
+ end if
+ prev_time = time
+!
+ if(time .ge. end_time) exit
+ end do
+ close(id_pick)
+!
+ acou = one / (end_time - true_start)
+ do ipick = 1, num_pick_sph_mode*num_pick_layer
+ do nd = 1, ntot_comp_pick_sph
+ ave_spec(nd,ipick) = ave_spec(nd,ipick) * acou
+ end do
+ end do
+!
+ call deallocate_pick_sph_monitor
+ call deallocate_num_pick_layer
+!
+! Evaluate standard deviation
+!
+ call open_sph_spec_read_monitor(id_pick)
+!
+ icou = 0
+ do
+ call read_sph_spec_4_monitor(id_pick, i_step, time, ierr)
+ if(ierr .gt. 0) exit
+!
+ if(time .ge. start_time) then
+ do ipick = 1, num_pick_sph_mode*num_pick_layer
+ do nd = 1, ntot_comp_pick_sph
+ prev_spec(nd,ipick) &
+ & = (d_rj_pick_sph_gl(nd,ipick) - ave_spec(nd,ipick))**2
+ end do
+ end do
+!
+ if(icou .eq. 0) then
+ true_start = time
+!
+ else
+ do ipick = 1, num_pick_sph_mode*num_pick_layer
+ do nd = 1, ntot_comp_pick_sph
+ sdev_spec(nd,ipick) = sdev_spec(nd,ipick) + half &
+ & * ((d_rj_pick_sph_gl(nd,ipick)-ave_spec(nd,ipick))**2 &
+ & + prev_spec(nd,ipick)) * (time - prev_time)
+ end do
+ end do
+ end if
+!
+ icou = icou + 1
+ write(*,*) 'step ', i_step, &
+ & ' is added for standard deviation: count is ', icou
+ end if
+ prev_time = time
+!
+ if(time .ge. end_time) exit
+ end do
+ close(id_pick)
+!
+ acou = one / (end_time - true_start)
+ do ipick = 1, num_pick_sph_mode*num_pick_layer
+ do nd = 1, ntot_comp_pick_sph
+ sdev_spec(nd,ipick) = sqrt(sdev_spec(nd,ipick)) * acou
+ end do
+ end do
+!
+! output time average
+!
+ do ipick = 1, num_pick_sph_mode*num_pick_layer
+ do nd = 1, ntot_comp_pick_sph
+ d_rj_pick_sph_gl(nd,ipick) = ave_spec(nd,ipick)
+ end do
+ end do
+!
+ pickup_sph_head = tave_header
+ call write_sph_spec_4_monitor(izero, i_step, time)
+!
+! output standard deviation
+!
+ do ipick = 1, num_pick_sph_mode*num_pick_layer
+ do nd = 1, ntot_comp_pick_sph
+ d_rj_pick_sph_gl(nd,ipick) = sdev_spec(nd,ipick)
+ end do
+ end do
+!
+ pickup_sph_head = sdev_header
+ call write_sph_spec_4_monitor(izero, i_step, time)
+!
+ call deallocate_pick_sph_monitor
+ call deallocate_num_pick_layer
+ deallocate(prev_spec, ave_spec, sdev_spec)
+!
+ write(*,*) '***** program finished *****'
+ stop
+!
+ end program tave_picked_sph_spec_data
More information about the CIG-COMMITS
mailing list