[cig-commits] [commit] Include_averaging_programs: Add Nusselt number averaging program (d367175)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon May 12 13:05:09 PDT 2014
Repository : https://github.com/geodynamics/calypso
On branch : Include_averaging_programs
Link : https://github.com/geodynamics/calypso/compare/0000000000000000000000000000000000000000...8e5b4fefdb5b27c6ed0b1448a44005ea47b02322
>---------------------------------------------------------------
commit d3671756c9a8dc30cd0c89470f2e2f0599bbe718
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Fri Apr 25 17:01:56 2014 -0700
Add Nusselt number averaging program
>---------------------------------------------------------------
d3671756c9a8dc30cd0c89470f2e2f0599bbe718
INSTALL | 3 +
README | 2 +-
doc/CALYPSO.pdf | Bin 2562754 -> 2567534 bytes
doc/CALYPSO.tex | 1 +
doc/tex_src/controls_CALYPSO.tex | 6 +
doc/tex_src/install_CALYPSO.tex | 3 +
doc/tex_src/outlines_CALYPSO.tex | 22 ++--
doc/tex_src/programs_CALYPSO.tex | 45 +++++++
.../Compositional_case_1/control_MHD | 2 +-
.../Compositional_case_1/control_snapshot | 2 +-
.../dynamobench_case_1/control_MHD | 10 +-
.../dynamobench_case_1/control_snapshot | 2 +-
.../dynamobench_case_2/control_MHD | 2 +-
.../dynamobench_case_2/control_snapshot | 2 +-
.../dynamo_benchmark/pseudo_vacuum/control_MHD | 2 +-
.../pseudo_vacuum/control_snapshot | 2 +-
examples/heat_composition_source/control_MHD | 2 +-
examples/heterogineous_temp/control_MHD | 2 +-
.../MHD_src/IO/sph_mhd_rms_IO.f90 | 7 +
.../MHD_src/sph_MHD/check_dependency_for_MHD.f90 | 1 +
.../MHD_src/sph_MHD/set_control_sph_mhd.f90 | 1 +
.../SPH_SHELL_src/pickup_gauss_coefficients.f90 | 57 +++++++-
.../SERIAL_src/IO/m_ctl_data_4_pickup_sph.f90 | 10 +-
.../SERIAL_src/IO/m_gauss_coefs_monitor_data.f90 | 8 +-
.../SERIAL_src/IO/m_no_heat_Nusselt_num.f90 | 144 +++++++++++++++++++++
.../SERIAL_src/IO/set_control_4_pickup_sph.f90 | 17 +++
src/programs/SPH_MHD/control_MHD | 2 +-
.../data_utilities/TIME_HISTORIES/Makefile | 20 ++-
.../TIME_HISTORIES/cal_tave_sph_ene_spectr.f90 | 16 +--
.../TIME_HISTORIES/m_sph_ene_spectra.f90 | 2 +-
.../TIME_HISTORIES/t_average_nusselt.f90 | 131 +++++++++++++++++++
.../TIME_HISTORIES/t_average_sph_ene_spec.f90 | 17 ++-
.../TIME_HISTORIES/tave_picked_sph_spec_data.f90 | 3 +-
33 files changed, 498 insertions(+), 48 deletions(-)
diff --git a/INSTALL b/INSTALL
index 25e8ae2..e92172b 100644
--- a/INSTALL
+++ b/INSTALL
@@ -2,6 +2,9 @@
! Brief installation procedure
!------------------------------------------------------------------------
!
+
+For detailed instllation procedure, please read documentation page 16--19.
+
The easiest way to install Calypso is the following:
1. Configuration
diff --git a/README b/README
index d6dd321..a7c3146 100644
--- a/README
+++ b/README
@@ -22,5 +22,5 @@ examples: Directory for examples
src: Directory for source files
work: Work directory (created by make)
-For installation instructions, please read the documentation in the manual or INSTALL.
+For installation instructions, please read INSTALL or the documentation page 16--19.
diff --git a/doc/CALYPSO.pdf b/doc/CALYPSO.pdf
index ebc362e..25edb08 100644
Binary files a/doc/CALYPSO.pdf and b/doc/CALYPSO.pdf differ
diff --git a/doc/CALYPSO.tex b/doc/CALYPSO.tex
index eb009ba..98bbc31 100644
--- a/doc/CALYPSO.tex
+++ b/doc/CALYPSO.tex
@@ -55,6 +55,7 @@ Calypso is a program package of magnetohydrodynamics (MHD) simulations in a rota
\newpage
\begin{appendices}
\input{tex_src/controls_CALYPSO.tex}
+\newpage
\input{tex_src/gpl3.tex}
\end{appendices}
diff --git a/doc/tex_src/controls_CALYPSO.tex b/doc/tex_src/controls_CALYPSO.tex
index 6553c4a..ab3be0e 100644
--- a/doc/tex_src/controls_CALYPSO.tex
+++ b/doc/tex_src/controls_CALYPSO.tex
@@ -617,6 +617,12 @@ File prefix for Gauss coefficients \verb|[gauss_coef_prefix]| is defined by Text
\verb|[gauss_coef_radius]| \\
Normalized radius to obtain Gauss coefficients \verb|[gauss_coef_radius]| is defined by real. Gauss coefficients are evaluated from the poloidal magnetic field at CMB by assuming electrically insulated mantle. Do not set \verb|[gauss_coef_radius]| less than the outer core radius $r_{o}$.
+\paragraph{\tt nusselt\_number\_prefix}
+\label{href_t:gauss_coefs_prefix}
+\verb|[nusselt_number_prefix]| \\
+File prefix for Nusselt number data at ICB and CMB \verb|[nusselt_number_prefix]| is defined by Text. Program add {\tt .dat} extension after this file prefix. If this file prefix is not defined, Nusselt number data are not generated. \\
+{\bf CAUTION: Nusselt number is not evaluated if heat source is exsist.}
+
\paragraph{\tt array pick\_layer\_ctl}
\label{href_t:pick_layer_ctl}
\verb|[Layer #]|
diff --git a/doc/tex_src/install_CALYPSO.tex b/doc/tex_src/install_CALYPSO.tex
index c652018..a38c2f0 100644
--- a/doc/tex_src/install_CALYPSO.tex
+++ b/doc/tex_src/install_CALYPSO.tex
@@ -144,6 +144,9 @@ Compile is performed using the {\tt make} command. The Makefile in the top direc
\item{\verb sph_dynamobench: } Data processing for dynamo benchmark test by Christensen {\it et. al.} (2002)
\item{\verb sph_zm_snapshot: } Generate zonal mean field
\item{\verb assemble_sph: } Data transfer program to change number of subdomains.
+\item{\verb t_ave_sph_mean_square: } Time averaging program for the mean square data.
+\item{\verb t_ave_picked_sph_coefs: } Time averaging program for the picked spectrum data.
+\item{\verb t_ave_nusselt: } Time averaging program for the Nusselt number data.
\item{\verb make_f90depends: } Program to generate dependency of the source code ({\tt make} command uses to generate {\tt work/Makefile})
\end{description}
%
diff --git a/doc/tex_src/outlines_CALYPSO.tex b/doc/tex_src/outlines_CALYPSO.tex
index d801106..c5b92e2 100644
--- a/doc/tex_src/outlines_CALYPSO.tex
+++ b/doc/tex_src/outlines_CALYPSO.tex
@@ -7,14 +7,17 @@ Calypso consists of programs shown in Table \ref{table:controls}.
\begin{tabular}{|c|c|c|}
\hline
Program & Control file name & Type \\ \hline
-\verb|gen_sph_grids| & \verb|control_sph_shell| & Serial \\
-\verb|sph_mhd| & \verb|control_MHD| & Parallel \\
-\verb|sph_initial_field| & \verb|control_MHD| & Parallel \\
-\verb|sph_add_initial_field| & \verb|control_MHD| & Parallel \\
-\verb|sph_snapshot| & \verb|control_snapshot| & Parallel \\
-\verb|sph_dynamobench| & \verb|control_snapshot| & Parallel \\
-\verb|sph_zm_snapshot| & \verb|control_snapshot| & Parallel \\
-\verb|assemble_sph| & \verb|control_sph_assemble| & Serial \\ \hline
+\verb|gen_sph_grids| & \verb|control_sph_shell| & Serial \\
+\verb|sph_mhd| & \verb|control_MHD| & Parallel \\ \hline
+\verb|sph_initial_field| & \verb|control_MHD| & Parallel \\
+\verb|sph_add_initial_field| & \verb|control_MHD| & Parallel \\
+\verb|sph_snapshot| & \verb|control_snapshot| & Parallel \\ \hline
+\verb|sph_zm_snapshot| & \verb|control_snapshot| & Parallel \\
+\verb|sph_dynamobench| & \verb|control_snapshot| & Parallel \\ \hline
+\verb|assemble_sph| & \verb|control_sph_assemble| & Serial \\ \hline
+\verb|t_ave_sph_mean_square| & N/A & Serial \\ \hline
+\verb|t_ave_picked_sph_coefs| & N/A & Serial \\ \hline
+\verb|t_ave_nusselt| & N/A & Serial \\ \hline
\end{tabular}
\end{center}
\label{table:controls}
@@ -46,9 +49,10 @@ To perform simulations by Calypso, the following processes are required.
\item Convert the parallel spectra data by \verb|assemble_sph| to continue with changing number of processes (if necessary).
\item Data analysis by \verb|sph_snapshot|, \verb|sph_snapshot|, or \verb|sph_dynamobench|.
\item Update initial fields by \verb|sph_add_initial_field| for more simulations (if necessary).
+\item Evaluate time averages by \verb|t_ave_sph_mean_square|, \verb|t_ave_picked_sph_coefs|, or \verb|t_ave_nusselt| if necessary.
\end{enumerate}
%
-The simulation program \verb|sph_mhd| requires an indexing file for spherical transform. \verb|sph_mhd| generates spectrum data and monitoring data, and field data in Cartesian coordinate as outputs. The data transform programs ( \verb|sph_snapshot| and \verb|sph_zm_snapshot|) generate outputs data from parallel spectra data. The flow of data is shown in Figure \ref{fig:flow_0}.
+The simulation program \verb|sph_mhd| requires an indexing file for spherical transform. \verb|sph_mhd| generates spectrum data and monitoring data, and field data in Cartesian coordinate as outputs. The data transform programs (\verb|sph_snapshot| and \verb|sph_zm_snapshot|) generate outputs data from parallel spectra data. The flow of data is shown in Figure \ref{fig:flow_0}.
%
\begin{figure}[H]
\begin{center}
diff --git a/doc/tex_src/programs_CALYPSO.tex b/doc/tex_src/programs_CALYPSO.tex
index ba22d57..7a93ad4 100644
--- a/doc/tex_src/programs_CALYPSO.tex
+++ b/doc/tex_src/programs_CALYPSO.tex
@@ -193,6 +193,7 @@ Data files for this program are listed in Table \ref{table:sph_mhd}. Indexing da
\verb|[layer_pwr_prefix]_lm.dat| & Single & Output \\ \hline
\verb|[gauss_coef_prefix].dat| & Single & Output \\
\verb|[picked_sph_prefix].dat| & Single & Output \\ \hline
+\verb|[nusselt_number_prefix].dat| & Single & Output \\ \hline
\verb|[fld_prefix].[step#].[domain#].[extension]| & - & Output \\ \hline
\end{tabular}
\end{center}
@@ -453,6 +454,8 @@ The format of the control file \verb|control_MHD| is described below. The detail
\hyperref[href_t:gauss_coefs_prefix]{(Detail)}
\item \verb|gauss_coefs_radius_ctl [gauss_coef_radius]| \\
\hyperref[href_t:gauss_coefs_radius_ctl]{(Detail)}
+ \item \verb|nusselt_number_prefix [nusselt_number_prefix]| \\
+ \hyperref[href_t:nusselt_number_prefix]{(Detail)}
\item \verb|array pick_layer_ctl [Layer #]|
\hyperref[href_t:pick_layer_ctl]{(Detail)}
\item \verb|array pick_sph_spectr_ctl [Degree] [Order]| \\
@@ -862,6 +865,22 @@ $
\end{description}
+\subsection{Nusselt number data {\tt [nusselt\_number\_prefix].dat}}
+{\bf CAUTION: Nusselt number is not evaluated if heat source is exsist.}
+The Nusselt number Nu at CMB and ICB is written for each step in one line. The Nusselt number is evaluated by
+%
+\begin{eqnarray*}
+Nu = \frac{<\partial T / \partial r>}{\partial T_{diff}/ \partial r},
+\end{eqnarray*}
+where, $<\partial T / \partial r>$ and $T_{diff}$ are the horizontal average of the temperature gradient at ICB and CMB and diffusive temperature profile, respectively. $T_{diff}$ is evaluated without heat source, as
+\begin{eqnarray*}
+T_{diff} = \frac{r_{o}T_{o} - r_{i}T_{i}}{r_{o} - r_{i}}
+ + \frac{r_{o}r_{i}\left(T_{i} - T_{o}\right)}{r_{o} - r_{i}} \frac{1}{r}.
+\end{eqnarray*}
+%
+This diffusive temperature profile is for the case without heat source in the fluid. If simulation is performed including the heat source, this data file does not written.
+If the Nusselt number data file exist before starting the simulation, programs append spectrum data at the end of files without checking constancy. If you change the configuration of data output structure, please move the old spectrum monitor file to another directory before starting the programs.
+
\newpage
\section{Data transform program \\
@@ -902,6 +921,7 @@ Simulation program outputs spectrum data as a whole field data. This program gen
\verb|[layer_pwr_prefix]_lm.dat| & Single & Output \\ \hline
\verb|[gauss_coef_prefix].dat| & Single & Output \\
\verb|[picked_sph_prefix].dat| & Single & Output \\ \hline
+\verb|[nusselt_number_prefix].dat| & Single & Output \\ \hline
\verb|[fld_prefix].[step#].[domain#].[extension]| & - & Output \\ \hline
\end{tabular}
\end{center}
@@ -1149,6 +1169,31 @@ This program is only used to generate Makefile in {\tt work} directory. Most of
\item The module name should be the same as the file name.
\end{itemize}
+\section{Time averaging programs ({\tt t\_ave\_picked\_sph\_coefs}, {\tt t\_ave\_sph\_mean\_square}, and {\tt t\_ave\_nusselt})}
+These small programs are used to evaluate time average and standard deviation of the time sequencial data.
+%
+\begin{table}[htdp]
+\caption{List of programs to take time average}
+\begin{center}
+\begin{tabular}{|c|c|}
+\hline
+ name & Program to use \\ \hline \hline
+\verb|[vol_pwr_prefix].dat| & \verb|t_ave_sph_mean_square| \\ \hline
+\verb|[vol_pwr_prefix]_l.dat| & \verb|t_ave_sph_mean_square| \\
+\verb|[vol_pwr_prefix]_m.dat| & \verb|t_ave_sph_mean_square| \\
+\verb|[vol_pwr_prefix]_lm.dat| & \verb|t_ave_sph_mean_square| \\ \hline
+\verb|[layer_pwr_prefix]_l.dat| & \verb|t_ave_sph_mean_square| \\
+\verb|[layer_pwr_prefix]_m.dat| & \verb|t_ave_sph_mean_square| \\
+\verb|[layer_pwr_prefix]_lm.dat| & \verb|t_ave_sph_mean_square| \\ \hline
+\verb|[picked_sph_prefix].dat| & \verb|t_ave_picked_sph_coefs| \\ \hline
+\verb|[nusselt_number_prefix].dat| & \verb|t_ave_sph_mean_square| \\ \hline
+\end{tabular}
+\end{center}
+\label{table:time_averages}
+\end{table}
+%
+
+
\section{Visualization using field data}
\label{sec:paraview}
The field data is written by XDMF or VTK data format using Cartesian coordinate. In this section we briefly introduce how to display the radial magnetic field using ParaView as an example.
diff --git a/examples/dynamo_benchmark/Compositional_case_1/control_MHD b/examples/dynamo_benchmark/Compositional_case_1/control_MHD
index 3848580..6a0b985 100644
--- a/examples/dynamo_benchmark/Compositional_case_1/control_MHD
+++ b/examples/dynamo_benchmark/Compositional_case_1/control_MHD
@@ -37,7 +37,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/examples/dynamo_benchmark/Compositional_case_1/control_snapshot b/examples/dynamo_benchmark/Compositional_case_1/control_snapshot
index cfade3c..00404c4 100644
--- a/examples/dynamo_benchmark/Compositional_case_1/control_snapshot
+++ b/examples/dynamo_benchmark/Compositional_case_1/control_snapshot
@@ -38,7 +38,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/examples/dynamo_benchmark/dynamobench_case_1/control_MHD b/examples/dynamo_benchmark/dynamobench_case_1/control_MHD
index 2de0c14..b466495 100644
--- a/examples/dynamo_benchmark/dynamobench_case_1/control_MHD
+++ b/examples/dynamo_benchmark/dynamobench_case_1/control_MHD
@@ -37,7 +37,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
@@ -237,8 +237,8 @@ begin MHD_control
!!!!! information for sime integration
begin time_step_ctl
elapsed_time_ctl 80000.
- i_step_init_ctl 10000
- i_step_finish_ctl 20000
+ i_step_init_ctl 0
+ i_step_finish_ctl 10000
!
i_step_check_ctl 10
i_step_rst_ctl 5000
@@ -250,8 +250,8 @@ begin MHD_control
!
!!!!! control for restart data
begin restart_file_ctl
-! rst_ctl dynamo_benchmark_1
- rst_ctl start_from_rst_file
+ rst_ctl dynamo_benchmark_1
+! rst_ctl start_from_rst_file
end restart_file_ctl
!
!!!!!! method for time evolution
diff --git a/examples/dynamo_benchmark/dynamobench_case_1/control_snapshot b/examples/dynamo_benchmark/dynamobench_case_1/control_snapshot
index 7735baf..3a91eea 100644
--- a/examples/dynamo_benchmark/dynamobench_case_1/control_snapshot
+++ b/examples/dynamo_benchmark/dynamobench_case_1/control_snapshot
@@ -38,7 +38,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/examples/dynamo_benchmark/dynamobench_case_2/control_MHD b/examples/dynamo_benchmark/dynamobench_case_2/control_MHD
index 33c25f6..745c0c4 100644
--- a/examples/dynamo_benchmark/dynamobench_case_2/control_MHD
+++ b/examples/dynamo_benchmark/dynamobench_case_2/control_MHD
@@ -37,7 +37,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/examples/dynamo_benchmark/dynamobench_case_2/control_snapshot b/examples/dynamo_benchmark/dynamobench_case_2/control_snapshot
index 8d84324..06f654d 100644
--- a/examples/dynamo_benchmark/dynamobench_case_2/control_snapshot
+++ b/examples/dynamo_benchmark/dynamobench_case_2/control_snapshot
@@ -37,7 +37,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/examples/dynamo_benchmark/pseudo_vacuum/control_MHD b/examples/dynamo_benchmark/pseudo_vacuum/control_MHD
index 7536711..55e5a8f 100644
--- a/examples/dynamo_benchmark/pseudo_vacuum/control_MHD
+++ b/examples/dynamo_benchmark/pseudo_vacuum/control_MHD
@@ -37,7 +37,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/examples/dynamo_benchmark/pseudo_vacuum/control_snapshot b/examples/dynamo_benchmark/pseudo_vacuum/control_snapshot
index 2a862a1..45449d9 100644
--- a/examples/dynamo_benchmark/pseudo_vacuum/control_snapshot
+++ b/examples/dynamo_benchmark/pseudo_vacuum/control_snapshot
@@ -37,7 +37,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/examples/heat_composition_source/control_MHD b/examples/heat_composition_source/control_MHD
index 9e705b1..988aa24 100644
--- a/examples/heat_composition_source/control_MHD
+++ b/examples/heat_composition_source/control_MHD
@@ -35,7 +35,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/examples/heterogineous_temp/control_MHD b/examples/heterogineous_temp/control_MHD
index dc88e4d..ddfbcc2 100644
--- a/examples/heterogineous_temp/control_MHD
+++ b/examples/heterogineous_temp/control_MHD
@@ -36,7 +36,7 @@ begin MHD_control
! vector_potential, magnetic_field, current_density, magnetic_potential
! composition, perturbation_temp
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/src/Fortran_libraries/MHD_src/IO/sph_mhd_rms_IO.f90 b/src/Fortran_libraries/MHD_src/IO/sph_mhd_rms_IO.f90
index 89ac3c5..3f87561 100644
--- a/src/Fortran_libraries/MHD_src/IO/sph_mhd_rms_IO.f90
+++ b/src/Fortran_libraries/MHD_src/IO/sph_mhd_rms_IO.f90
@@ -47,8 +47,10 @@
subroutine output_rms_sph_mhd_control
!
use m_t_step_parameter
+ use m_boundary_params_sph_MHD
use set_exit_flag_4_visualizer
use cal_rms_fields_by_sph
+ use m_no_heat_Nusselt_num
!
integer (kind = kint) :: i_flag
!
@@ -60,6 +62,8 @@
call cal_rms_sph_outer_core
call cal_gauss_coefficients
call pickup_sph_spec_4_monitor
+ call cal_no_heat_source_Nu(sph_bc_U%kr_in, sph_bc_U%kr_out, &
+ & sph_bc_U%r_ICB(0), sph_bc_U%r_CMB(0) )
!
call write_total_energy_to_screen(my_rank, i_step_MHD, time)
!
@@ -71,6 +75,9 @@
call write_gauss_coefs_4_monitor(my_rank, istep_max_dt, time)
call write_sph_spec_4_monitor(my_rank, istep_max_dt, time)
!
+ call write_no_heat_source_Nu(idx_rj_degree_zero, &
+ & istep_max_dt, time)
+!
end subroutine output_rms_sph_mhd_control
!
! --------------------------------------------------------------------
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/check_dependency_for_MHD.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/check_dependency_for_MHD.f90
index c7f610b..92f9f37 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/check_dependency_for_MHD.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/check_dependency_for_MHD.f90
@@ -97,6 +97,7 @@
& .or. phys_nod_name(i) .eq. fhd_filter_temp &
& .or. phys_nod_name(i) .eq. fhd_buoyancy &
& .or. phys_nod_name(i) .eq. fhd_heat_source &
+ & .or. phys_nod_name(i) .eq. fhd_grad_temp &
& ) then
num_check = 1
phys_check_name(1) = fhd_temp
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_control_sph_mhd.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_control_sph_mhd.f90
index 03406ad..118af6f 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_control_sph_mhd.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_control_sph_mhd.f90
@@ -120,6 +120,7 @@
!
call set_ctl_params_pick_sph
call set_ctl_params_pick_gauss
+ call set_ctl_params_no_heat_Nu
!
end subroutine set_control_4_sph_mhd
!
diff --git a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_gauss_coefficients.f90 b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_gauss_coefficients.f90
index b68067f..e3f68d1 100644
--- a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_gauss_coefficients.f90
+++ b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/pickup_gauss_coefficients.f90
@@ -5,10 +5,13 @@
!!@date Programmed in Dec., 2012
!
!> @brief choose Gauss coefficients to output
+!>@n Evaluate Nusselt number without heat source
!!
!!@verbatim
!! subroutine init_gauss_coefs_4_monitor
!! subroutine cal_gauss_coefficients
+!!
+!! subroutine cal_no_heat_source_Nu(kr_ICB, kr_CMB, r_in, r_out)
!!@endverbatim
!
module pickup_gauss_coefficients
@@ -17,7 +20,6 @@
use m_constants
!
use m_spheric_parameter
- use m_gauss_coefs_monitor_data
use pickup_sph_spectr
!
implicit none
@@ -33,6 +35,7 @@
subroutine init_gauss_coefs_4_monitor
!
use m_sph_phys_address
+ use m_gauss_coefs_monitor_data
!
!
if (ipol%i_magne .gt. 0) then
@@ -71,6 +74,7 @@
use calypso_mpi
use m_sph_spectr_data
use m_sph_phys_address
+ use m_gauss_coefs_monitor_data
!
integer(kind = kint) :: inum, j, l, inod
real(kind = kreal) :: a2r_4_gauss
@@ -110,6 +114,7 @@
subroutine set_gauss_coefs_labels
!
use set_parallel_file_name
+ use m_gauss_coefs_monitor_data
!
integer(kind = kint) :: j, l, m, mm, inum
character(len=kchara) :: gauss_head
@@ -135,5 +140,55 @@
end subroutine set_gauss_coefs_labels
!
! -----------------------------------------------------------------------
+! -----------------------------------------------------------------------
+!
+ subroutine cal_no_heat_source_Nu(kr_ICB, kr_CMB, r_in, r_out)
+!
+ use m_sph_spectr_data
+ use m_sph_phys_address
+ use m_no_heat_Nusselt_num
+!
+ integer(kind = kint), intent(in) :: kr_ICB, kr_CMB
+ real(kind = kreal), intent(in) :: r_in, r_out
+ real(kind = kreal) :: temp_ICB, dTdr_ICB
+ real(kind = kreal) :: temp_CMB, dTdr_CMB
+!
+ real(kind = kreal) :: c1, c2
+! real(kind = kreal) :: dTdr_diff_ICB, dTdr_diff_CMB
+ integer(kind = kint) :: inod_ICB, inod_CMB
+!
+!
+ if(iflag_no_source_Nu .eq. izero) return
+ if(idx_rj_degree_zero .eq. 0) return
+!
+ r_ICB_Nu = r_in
+ r_CMB_Nu = r_out
+!
+ inod_ICB = idx_rj_degree_zero + (kr_ICB-1) * nidx_rj(2)
+ temp_ICB = d_rj(inod_ICB,ipol%i_temp)
+! dTdr_ICB = half*d_rj(inod_ICB,ipol%i_grad_t) &
+! & * a_r_1d_rj_r(kr_ICB)**2
+!
+ inod_CMB = idx_rj_degree_zero + (kr_CMB-1) * nidx_rj(2)
+ temp_CMB = d_rj(inod_CMB,ipol%i_temp)
+! dTdr_CMB = half*d_rj(inod_CMB,ipol%i_grad_t) &
+! & * a_r_1d_rj_r(kr_CMB)**2
+!
+ c1 = (r_CMB_Nu*temp_CMB - r_ICB_Nu*temp_ICB) &
+ & / ( r_CMB_Nu - r_ICB_Nu )
+ c2 = r_CMB_Nu * r_ICB_Nu * (temp_ICB - temp_CMB) &
+ & / ( r_CMB_Nu - r_ICB_Nu )
+!
+! dTdr_diff_ICB = - c2 * a_r_1d_rj_r(kr_ICB)**2
+! dTdr_diff_CMB = - c2 * a_r_1d_rj_r(kr_CMB)**2
+! Nu_ICB = dTdr_ICB / dTdr_diff_ICB
+! Nu_CMB = dTdr_CMB / dTdr_diff_CMB
+!
+ Nu_ICB = - half*d_rj(inod_ICB,ipol%i_grad_t) / c2
+ Nu_CMB = - half*d_rj(inod_CMB,ipol%i_grad_t) / c2
+!
+ end subroutine cal_no_heat_source_Nu
+!
+! -----------------------------------------------------------------------
!
end module pickup_gauss_coefficients
diff --git a/src/Fortran_libraries/SERIAL_src/IO/m_ctl_data_4_pickup_sph.f90 b/src/Fortran_libraries/SERIAL_src/IO/m_ctl_data_4_pickup_sph.f90
index 330ff2a..a044e32 100644
--- a/src/Fortran_libraries/SERIAL_src/IO/m_ctl_data_4_pickup_sph.f90
+++ b/src/Fortran_libraries/SERIAL_src/IO/m_ctl_data_4_pickup_sph.f90
@@ -31,6 +31,7 @@
!!
!! picked_sph_prefix 'sph_spectr/picked_mode'
!! gauss_coefs_prefix 'sph_spectr/gauss_coefs'
+!! nusselt_number_prefix 'Nusselt'
!!
!! gauss_coefs_radius_ctl 2.91
!!
@@ -94,6 +95,7 @@
character(len = kchara) :: volume_pwr_spectr_prefix
character(len = kchara) :: layered_pwr_spectr_prefix
!
+ character(len = kchara) :: Nusselt_file_prefix
character(len = kchara) :: picked_mode_head_ctl
character(len = kchara) :: gauss_coefs_prefix
real(kind = kreal) :: gauss_coefs_radius_ctl = 2.91
@@ -144,6 +146,8 @@
character(len=kchara), parameter &
& :: hd_gauss_coefs_head = 'gauss_coefs_prefix'
character(len=kchara), parameter &
+ & :: hd_Nusselt_file_head = 'nusselt_number_prefix'
+ character(len=kchara), parameter &
& :: hd_gauss_coefs_r = 'gauss_coefs_radius_ctl'
!
character(len=kchara), parameter &
@@ -178,6 +182,7 @@
integer (kind=kint) :: i_layer_rms_head = 0
integer (kind=kint) :: i_picked_mode_head = 0
integer (kind=kint) :: i_gauss_coefs_head = 0
+ integer (kind=kint) :: i_Nusselt_file_head = 0
integer (kind=kint) :: i_gauss_coefs_r = 0
!
integer (kind=kint) :: i_num_pick_layer = 0
@@ -197,7 +202,7 @@
!
private :: hd_pick_sph, i_pick_sph, hd_num_pick_layer
private :: hd_gauss_coefs_head, hd_gauss_coefs_r
- private :: hd_picked_mode_head
+ private :: hd_picked_mode_head, hd_Nusselt_file_head
private :: hd_num_pick_sph, hd_num_pick_gauss
private :: hd_num_pick_l, hd_num_pick_m
private :: hd_num_pick_gauss_l, hd_num_pick_gauss_m
@@ -426,6 +431,9 @@
call read_character_ctl_item(hd_picked_mode_head, &
& i_picked_mode_head, picked_mode_head_ctl)
!
+ call read_character_ctl_item(hd_Nusselt_file_head, &
+ & i_Nusselt_file_head, Nusselt_file_prefix)
+!
call read_character_ctl_item(hd_voume_ave_head, &
& i_voume_ave_head, volume_average_prefix)
call read_character_ctl_item(hd_voume_rms_head, &
diff --git a/src/Fortran_libraries/SERIAL_src/IO/m_gauss_coefs_monitor_data.f90 b/src/Fortran_libraries/SERIAL_src/IO/m_gauss_coefs_monitor_data.f90
index 043c470..82be2ae 100644
--- a/src/Fortran_libraries/SERIAL_src/IO/m_gauss_coefs_monitor_data.f90
+++ b/src/Fortran_libraries/SERIAL_src/IO/m_gauss_coefs_monitor_data.f90
@@ -40,8 +40,6 @@
integer(kind = kint), parameter :: id_gauss_coef = 23
!> File prefix for Gauss coefficients file
character(len = kchara) :: gauss_coefs_file_head
-!> File name for Gauss coefficients file
- character(len = kchara) :: gauss_coefs_file_name
!
!> Number of modes of Gauss coefficients to be evaluated
integer(kind = kint) :: num_pick_gauss_coefs = 0
@@ -158,6 +156,8 @@
use set_parallel_file_name
use write_field_labels
!
+ character(len = kchara) :: gauss_coefs_file_name
+!
!
call add_dat_extension(gauss_coefs_file_head, &
& gauss_coefs_file_name)
@@ -216,13 +216,17 @@
subroutine open_gauss_coefs_read_monitor(id_pick)
!
use skip_comment_f
+ use set_parallel_file_name
!
integer(kind = kint), intent(in) :: id_pick
!
integer(kind = kint) :: i
+ character(len = kchara) :: gauss_coefs_file_name
character(len=255) :: tmpchara
!
!
+ call add_dat_extension(gauss_coefs_file_head, &
+ & gauss_coefs_file_name)
open(id_pick, file = gauss_coefs_file_name)
!
call skip_comment(tmpchara,id_pick)
diff --git a/src/Fortran_libraries/SERIAL_src/IO/m_no_heat_Nusselt_num.f90 b/src/Fortran_libraries/SERIAL_src/IO/m_no_heat_Nusselt_num.f90
new file mode 100644
index 0000000..289838c
--- /dev/null
+++ b/src/Fortran_libraries/SERIAL_src/IO/m_no_heat_Nusselt_num.f90
@@ -0,0 +1,144 @@
+!>@file m_no_heat_Nusselt_num.f90
+!!@brief module m_no_heat_Nusselt_num
+!!
+!!@author H. Matsui
+!!@date Programmed in Dec., 2012
+!
+!>@brief Data arrays for Nusselt number
+!!
+!!@verbatim
+!! subroutine write_no_heat_source_Nu(idx_rj_degree_zero, &
+!! & i_step, time)
+!!
+!! subroutine open_read_no_heat_source_Nu(id_pick)
+!! subroutine read_no_heat_source_Nu(id_pick, i_step, time, ierr)
+!!@endverbatim
+!
+ module m_no_heat_Nusselt_num
+!
+ use m_precision
+ use m_constants
+!
+ implicit none
+!
+!
+!> Output flag for Nusselt number IO
+ integer(kind = kint) :: iflag_no_source_Nu = 0
+!> File ID for Nusselt number IO
+ integer(kind = kint), parameter :: id_Nusselt = 23
+!> File prefix for Nusselt number file
+ character(len = kchara) :: Nusselt_file_head = 'Nusselt'
+!
+!> Radius at inner boundary
+ real(kind = kreal) :: r_ICB_Nu
+!> Radius at outer boundary
+ real(kind = kreal) :: r_CMB_Nu
+!> Nusselt number at inner boundary
+ real(kind = kreal) :: Nu_ICB
+!> Nusselt number at outer boundary
+ real(kind = kreal) :: Nu_CMB
+!
+ private :: open_no_heat_source_Nu
+!
+! -----------------------------------------------------------------------
+!
+ contains
+!
+! -----------------------------------------------------------------------
+!
+ subroutine open_no_heat_source_Nu
+!
+ use set_parallel_file_name
+ use write_field_labels
+!
+ character(len = kchara) :: Nusselt_file_name
+!
+!
+ call add_dat_extension(Nusselt_file_head, Nusselt_file_name)
+ open(id_Nusselt, file = Nusselt_file_name, &
+ & form='formatted', status='old', position='append', err = 99)
+ return
+!
+ 99 continue
+ open(id_Nusselt, file = Nusselt_file_name, &
+ & form='formatted', status='replace')
+!
+!
+ write(id_Nusselt,'(a)') '# Inner_radius, Outer_radius'
+ write(id_Nusselt,'(1p2e25.15e3)') r_ICB_Nu, r_CMB_Nu
+!
+ write(id_Nusselt,'(a)',advance='NO') &
+ & 't_step time Nu_ICB Nu_CMB'
+ write(id_Nusselt,'(a)') ''
+!
+ end subroutine open_no_heat_source_Nu
+!
+! -----------------------------------------------------------------------
+!
+ subroutine write_no_heat_source_Nu(idx_rj_degree_zero, &
+ & i_step, time)
+!
+ integer(kind = kint), intent(in) :: idx_rj_degree_zero, i_step
+ real(kind = kreal), intent(in) :: time
+!
+!
+ if(iflag_no_source_Nu .eq. izero) return
+ if(idx_rj_degree_zero .eq. izero) return
+!
+ call open_no_heat_source_Nu
+!
+ write(id_Nusselt,'(i10,1p3e23.14e3)') &
+ & i_step, time, Nu_ICB, Nu_CMB
+!
+ close(id_Nusselt)
+!
+ end subroutine write_no_heat_source_Nu
+!
+! -----------------------------------------------------------------------
+! -----------------------------------------------------------------------
+!
+ subroutine open_read_no_heat_source_Nu(id_pick)
+!
+ use skip_comment_f
+ use set_parallel_file_name
+!
+ integer(kind = kint), intent(in) :: id_pick
+ character(len = kchara) :: Nusselt_file_name
+!
+ integer(kind = kint) :: i
+ character(len=255) :: tmpchara
+!
+!
+ call add_dat_extension(Nusselt_file_head, Nusselt_file_name)
+ open(id_pick, file = Nusselt_file_name)
+!
+ call skip_comment(tmpchara,id_pick)
+ read(tmpchara,*) r_ICB_Nu, r_CMB_Nu
+!
+ read(id_pick,*) (tmpchara,i=1,4)
+!
+ end subroutine open_read_no_heat_source_Nu
+!
+! -----------------------------------------------------------------------
+!
+ subroutine read_no_heat_source_Nu(id_pick, i_step, time, ierr)
+!
+ integer(kind = kint), intent(in) :: id_pick
+ integer(kind = kint), intent(inout) :: i_step, ierr
+ real(kind = kreal), intent(inout) :: time
+!
+!
+ ierr = 0
+ read(id_pick,*,err=99,end=99) i_step, time, Nu_ICB, Nu_CMB
+!
+ return
+!
+ 99 continue
+ ierr = 1
+ return
+!
+ end subroutine read_no_heat_source_Nu
+!
+! -----------------------------------------------------------------------
+!
+ end module m_no_heat_Nusselt_num
diff --git a/src/Fortran_libraries/SERIAL_src/IO/set_control_4_pickup_sph.f90 b/src/Fortran_libraries/SERIAL_src/IO/set_control_4_pickup_sph.f90
index cc13855..43a8006 100644
--- a/src/Fortran_libraries/SERIAL_src/IO/set_control_4_pickup_sph.f90
+++ b/src/Fortran_libraries/SERIAL_src/IO/set_control_4_pickup_sph.f90
@@ -171,4 +171,21 @@
!
! -----------------------------------------------------------------------
!
+ subroutine set_ctl_params_no_heat_Nu
+!
+ use m_ctl_data_4_pickup_sph
+ use m_no_heat_Nusselt_num
+!
+!
+ if(i_Nusselt_file_head .gt. 0) then
+ iflag_no_source_Nu = 1
+ Nusselt_file_head = Nusselt_file_prefix
+ else
+ iflag_no_source_Nu = 0
+ end if
+!
+ end subroutine set_ctl_params_no_heat_Nu
+!
+! -----------------------------------------------------------------------
+!
end module set_control_4_pickup_sph
diff --git a/src/programs/SPH_MHD/control_MHD b/src/programs/SPH_MHD/control_MHD
index a558e59..d18823c 100644
--- a/src/programs/SPH_MHD/control_MHD
+++ b/src/programs/SPH_MHD/control_MHD
@@ -49,7 +49,7 @@ begin MHD_control
! kinetic_helicity, magnetic_helicity
! current_helicity, cross_helicity
!
-! buoyancy_flux, Lorentz_work, mag_tension_work
+! buoyancy_flux, Lorentz_work
! magnetic_ene_generation, work_against_Lorentz
! temp_generation, part_temp_gen
! vis_ene_diffuse, mag_ene_diffuse
diff --git a/src/programs/data_utilities/TIME_HISTORIES/Makefile b/src/programs/data_utilities/TIME_HISTORIES/Makefile
index 259d683..445ac1e 100644
--- a/src/programs/data_utilities/TIME_HISTORIES/Makefile
+++ b/src/programs/data_utilities/TIME_HISTORIES/Makefile
@@ -6,6 +6,7 @@ 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
+TARGET_T_AVE_NUSSELT = t_ave_nusselt
SOURCES = $(shell ls *.f90)
@@ -18,6 +19,9 @@ cal_tave_sph_ene_spectr.o
MOD_T_AVE_PICK_SPH = \
tave_picked_sph_spec_data.o
+MOD_T_AVE_NUSSELT = \
+t_average_nusselt.o
+
#
# ------------------------------------------------------------------
#
@@ -30,12 +34,14 @@ target_list:
>> $(MAKENAME)
@echo 'TARGET_T_AVE_PICK_SPH_SPEC = $$(BUILDDIR)/$(TARGET_T_AVE_PICK_SPH_SPEC)' \
>> $(MAKENAME)
+ @echo 'TARGET_T_AVE_NUSSELT = $$(BUILDDIR)/$(TARGET_T_AVE_NUSSELT)' \
+ >> $(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)
+ '$$(TARGET_TAVE_SPH_ENE_SPEC) $$(TARGET_T_AVE_NUSSELT)' \
+ '$$(TARGET_T_AVE_PICK_SPH_SPEC) ' >> $(MAKENAME)
@echo >> $(MAKENAME)
@echo '$$(TARGET_TAVE_SPH_ENE_SPEC)': '$$(MOD_TAVE_SPH_ENE_SPEC)' \
'$$(LIB_FILES_SPH_MHD)' \
@@ -53,6 +59,14 @@ target_task:
'-L. $$(LIBS_SPH_MHD)' \
'$$(F90LIBS)' >> $(MAKENAME)
@echo '' >> $(MAKENAME)
+ @echo '$$(TARGET_T_AVE_NUSSELT)': '$$(MOD_T_AVE_NUSSELT)' \
+ '$$(LIB_FILES_SPH_MHD)' \
+ >> $(MAKENAME)
+ @echo ' ''$$(F90)' '$$(F90FLAGS)' -o '$$(TARGET_T_AVE_NUSSELT)' \
+ '$$(MOD_T_AVE_NUSSELT)' \
+ '-L. $$(LIBS_SPH_MHD)' \
+ '$$(F90LIBS)' >> $(MAKENAME)
+ @echo '' >> $(MAKENAME)
lib_name:
@@ -62,6 +76,8 @@ mod_list:
@echo $(MOD_TAVE_SPH_ENE_SPEC) >> $(MAKENAME)
@echo MOD_T_AVE_PICK_SPH= \\ >> $(MAKENAME)
@echo $(MOD_T_AVE_PICK_SPH) >> $(MAKENAME)
+ @echo MOD_T_AVE_NUSSELT= \\ >> $(MAKENAME)
+ @echo $(MOD_T_AVE_NUSSELT) >> $(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
index 85f207e..3c712ea 100644
--- 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
@@ -53,7 +53,7 @@
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)
+ & * (time_sph - pre_time)
spectr_pre_t(nd,kr) = tmp
end do
!$omp end do nowait
@@ -67,13 +67,13 @@
!
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)
+ & * (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)
+ & * (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)
+ & * (time_sph - pre_time)
!
spectr_pre_l(nd,lth,kr) = tmp_l
spectr_pre_m(nd,lth,kr) = tmp_m
@@ -103,7 +103,7 @@
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)
+ & * (time_sph - pre_time)
spectr_pre_t(nd,kr) = tmp
end do
!$omp end do nowait
@@ -117,13 +117,13 @@
!
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)
+ & * (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)
+ & * (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)
+ & * (time_sph - pre_time)
!
spectr_pre_l(nd,lth,kr) = tmp_l
spectr_pre_m(nd,lth,kr) = tmp_m
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
index c84354a..221d35c 100644
--- a/src/programs/data_utilities/TIME_HISTORIES/m_sph_ene_spectra.f90
+++ b/src/programs/data_utilities/TIME_HISTORIES/m_sph_ene_spectra.f90
@@ -134,7 +134,7 @@
subroutine select_sph_ene_spec_data_file
!
!
- write(*,*) ' Choose data type to taking average'
+ write(*,*) ' Choose data type to take average'
write(*,*) ' 1: volume average spectrum '
write(*,*) ' 2: spectrum on each layer '
write(*,*) ' 3: spectrum on specific layer '
diff --git a/src/programs/data_utilities/TIME_HISTORIES/t_average_nusselt.f90 b/src/programs/data_utilities/TIME_HISTORIES/t_average_nusselt.f90
new file mode 100644
index 0000000..d795dcc
--- /dev/null
+++ b/src/programs/data_utilities/TIME_HISTORIES/t_average_nusselt.f90
@@ -0,0 +1,131 @@
+!t_average_nusselt.f90
+! program t_average_nusselt
+!
+! programmed by H.Matsui on Apr., 2014
+!
+ program t_average_nusselt
+!
+ use m_precision
+ use m_constants
+!
+ use m_no_heat_Nusselt_num
+!
+ implicit none
+!
+ real(kind = kreal) :: prev_Nu(2)
+ real(kind = kreal) :: ave_Nu(2)
+ real(kind = kreal) :: sdev_Nu(2)
+ integer(kind = kint), parameter :: id_pick = 15
+!
+ integer(kind = kint) :: i_step, ierr, icou, ipick, nd, i
+ integer(kind = kint) :: istep_start
+ real(kind = kreal) :: acou, time, prev_time
+ real(kind = kreal) :: start_time, end_time, true_start
+!
+!
+ write(*,*) '-----------------------------------------------'
+ write(*,*) ' CAUTION!!'
+ write(*,*) 'This program can evaluate Nusselt number'
+ write(*,*) 'only when there is NO internal heat source'
+ write(*,*) 'in the outer core!!'
+ write(*,*) '-----------------------------------------------'
+ write(*,*) ''
+ write(*,*) 'Input picked harmonics coefficients file prefix'
+ read(5,*) Nusselt_file_head
+!
+ write(*,*) 'Input start and end time'
+ read(5,*) start_time, end_time
+!
+ call open_read_no_heat_source_Nu(id_pick)
+!
+! Evaluate time average
+!
+ icou = 0
+ write(*,'(a5,i12,a30,i12)',advance="NO") &
+ & 'step= ', i_step, ' averaging finished. Count= ', icou
+ do
+ call read_no_heat_source_Nu(id_pick, i_step, time, ierr)
+ if(ierr .gt. 0) exit
+!
+ if(time .ge. start_time) then
+ prev_Nu(1) = Nu_ICB
+ prev_Nu(2) = Nu_CMB
+!
+ if(icou .eq. 0) then
+ true_start = time
+ else
+ ave_Nu(1) = ave_Nu(1) &
+ & + half*(Nu_ICB + prev_Nu(1)) * (time - prev_time)
+ ave_Nu(2) = ave_Nu(2) &
+ & + half*(Nu_CMB + prev_Nu(2)) * (time - prev_time)
+ end if
+!
+ icou = icou + 1
+ end if
+ prev_time = time
+!
+ write(*,'(59a1,a5,i12,a30,i12)',advance="NO") (char(8),i=1,59), &
+ & 'step= ', i_step, ' averaging finished. Count= ', icou
+ if(time .ge. end_time) exit
+ end do
+ write(*,*)
+ close(id_pick)
+!
+ acou = one / (end_time - true_start)
+ ave_Nu(1:2) = ave_Nu(1:2) * acou
+!
+! Evaluate standard deviation
+!
+ call open_read_no_heat_source_Nu(id_pick)
+!
+ write(*,'(a5,i12,a30,i12)',advance="NO") &
+ & 'step= ', i_step, ' deviation finished. Count= ', icou
+ icou = 0
+ do
+ call read_no_heat_source_Nu(id_pick, i_step, time, ierr)
+ if(ierr .gt. 0) exit
+!
+ if(time .ge. start_time) then
+ prev_Nu(1) = (Nu_ICB - ave_Nu(1))**2
+ prev_Nu(2) = (Nu_CMB - ave_Nu(2))**2
+!
+ if(icou .eq. 0) then
+ true_start = time
+ else
+ sdev_Nu(1) = sdev_Nu(1) &
+ & + half*( (Nu_ICB - ave_Nu(1))**2 + prev_Nu(1)) &
+ & * (time - prev_time)
+ sdev_Nu(2) = sdev_Nu(2) &
+ & + half*( (Nu_CMB - ave_Nu(2))**2 + prev_Nu(2)) &
+ & * (time - prev_time)
+ end if
+!
+ icou = icou + 1
+ end if
+ prev_time = time
+!
+ write(*,'(59a1,a5,i12,a30,i12)',advance="NO") (char(8),i=1,59), &
+ & 'step= ', i_step, ' deviation finished. Count= ', icou
+ if(time .ge. end_time) exit
+ end do
+ write(*,*)
+ close(id_pick)
+!
+ acou = one / (end_time - true_start)
+ sdev_Nu(1:2) = sqrt(sdev_Nu(1:2)) * acou
+!
+! output Results
+!
+ write(*,'(a,1p2e25.15e3)') 'Inner and outer radius: ', &
+ & r_ICB_Nu, r_CMB_Nu
+ write(*,'(a,1p2e25.15e3)') 'Start and end time: ', &
+ & true_start, end_time
+ write(*,'(a)') 'Average and Std. Dev. of Nu at ICB:'
+ write(*,'(1p2e25.15e3)') ave_Nu(1), sdev_Nu(1)
+ write(*,'(a)') 'Average and Std. Dev. of Nu at CMB:'
+ write(*,'(1p2e25.15e3)') ave_Nu(2), sdev_Nu(2)
+!
+ write(*,*) '***** program finished *****'
+ stop
+!
+ end program t_average_nusselt
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
index bc7694b..bdc803c 100644
--- 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
@@ -18,7 +18,7 @@
implicit none
!
!
- integer(kind = kint) :: ierr
+ integer(kind = kint) :: ierr, i
integer(kind = kint) :: icou, istep
real(kind = kreal) :: start_time, end_time
!
@@ -45,6 +45,8 @@
!
ist_true = -1
icou = 0
+ write(*,'(a5,i12,a30,i12)',advance="NO") &
+ & 'step= ', istep, ' averaging finished. Count= ', icou
do
if(iflag_sph_ene_file .eq. 1) then
call read_org_volume_ene_data(istep, ierr)
@@ -67,12 +69,13 @@
end if
end if
!
+ write(*,'(59a1,a5,i12,a30,i12)',advance="NO") (char(8),i=1,59), &
+ & 'step= ', istep, ' averaging finished. Count= ', icou
if (time_sph .ge. end_time) exit
-!
- write(*,*) 'step', istep, 'averaging finished. Count: ', icou
end do
99 continue
call close_ene_spec_data
+ write(*,*)
!
call divide_average_ene_sph
call output_tave_ene_sph_data
@@ -85,6 +88,8 @@
!
ist_true = -1
icou = 0
+ write(*,'(a5,i12,a30,i12)',advance="NO") &
+ & 'step= ', istep, ' deviation finished. Count= ', icou
do
if(iflag_sph_ene_file .eq. 1) then
call read_org_volume_ene_data(istep, ierr)
@@ -107,12 +112,12 @@
end if
end if
!
+ write(*,'(59a1,a5,i12,a30,i12)',advance="NO") (char(8),i=1,59), &
+ & 'step= ', istep, ' deviation finished. Count= ', icou
if (time_sph .ge. end_time) exit
-!
- write(*,*) 'step', istep, 'deviation finished. Count: ', icou
-!
end do
98 continue
+ write(*,*)
call close_ene_spec_data
!
call divide_deviation_ene_sph
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
index f28d0d9..aaa930b 100644
--- 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
@@ -21,9 +21,8 @@
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
+ real(kind = kreal) :: start_time, end_time, true_start
!
!
write(*,*) 'Input picked spectr evolution file header'
More information about the CIG-COMMITS
mailing list