[cig-commits] [commit] master: Write out all model parameters into NetCDF file (dd6f812)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Sat Sep 6 07:02:51 PDT 2014
Repository : https://github.com/geodynamics/axisem
On branch : master
Link : https://github.com/geodynamics/axisem/compare/b2c8f921d7cb70dd1d0edb25c35f6bd6de572cc7...dd6f81229fe57323788baa40ff8edb224d1e10eb
>---------------------------------------------------------------
commit dd6f81229fe57323788baa40ff8edb224d1e10eb
Author: Simon Stähler <staehler at geophysik.uni-muenchen.de>
Date: Sat Sep 6 15:56:43 2014 +0200
Write out all model parameters into NetCDF file
- Now includes: phi, xi, eta and Qmu and Qkappa
- Fix in nc_dump_elastic_parameters, where solid elements were written
into fluid part.
>---------------------------------------------------------------
dd6f81229fe57323788baa40ff8edb224d1e10eb
SOLVER/get_model.F90 | 21 ++++++--
SOLVER/nc_routines.F90 | 135 ++++++++++++++++++++++++++++++++++++++++++++++---
2 files changed, 143 insertions(+), 13 deletions(-)
diff --git a/SOLVER/get_model.F90 b/SOLVER/get_model.F90
index 617a928..692baaf 100644
--- a/SOLVER/get_model.F90
+++ b/SOLVER/get_model.F90
@@ -79,7 +79,7 @@ module get_model
!! xmgrace timestep_rad.dat or xmgrace period_rad.dat
!-----------------------------------------------------------------------------
subroutine read_model(rho, lambda, mu, xi_ani, phi_ani, eta_ani, &
- fa_ani_theta, fa_ani_phi, Q_mu, Q_kappa)
+ fa_ani_theta, fa_ani_phi, Q_mu_1d, Q_kappa_1d)
use commun, ONLY : barrier
use lateral_heterogeneities
@@ -92,8 +92,9 @@ subroutine read_model(rho, lambda, mu, xi_ani, phi_ani, eta_ani, &
real(kind=dp), dimension(0:npol,0:npol,nelem), intent(out) :: xi_ani, phi_ani, eta_ani
real(kind=dp), dimension(0:npol,0:npol,nelem), intent(out) :: fa_ani_theta, fa_ani_phi
- real(kind=realkind), dimension(nel_solid), intent(out), optional :: Q_mu, Q_kappa
+ real(kind=realkind), dimension(nel_solid), intent(out), optional :: Q_mu_1d, Q_kappa_1d
+ real(kind=dp), dimension(0:npol,0:npol,nelem) :: Q_mu, Q_kappa
real(kind=dp) :: s,z,r,theta,r1
real(kind=dp) :: vphtmp, vpvtmp, vshtmp, vsvtmp
integer :: iel,ipol,jpol,iidom,ieldom(nelem),domcount(ndisc),iel_count
@@ -244,8 +245,8 @@ subroutine read_model(rho, lambda, mu, xi_ani, phi_ani, eta_ani, &
do iel=1, nel_solid
iidom = ieldom(ielsolid(iel))
call compute_coordinates(s, z, r, theta, ielsolid(iel), npol/2 - 1, npol/2 - 1)
- Q_mu(iel) = velocity(r, 'Qmu', iidom, modelstring, lfbkgrdmodel)
- Q_kappa(iel) = velocity(r, 'Qka', iidom, modelstring, lfbkgrdmodel)
+ Q_mu_1d(iel) = velocity(r, 'Qmu', iidom, modelstring, lfbkgrdmodel)
+ Q_kappa_1d(iel) = velocity(r, 'Qka', iidom, modelstring, lfbkgrdmodel)
enddo
endif
@@ -263,7 +264,7 @@ subroutine read_model(rho, lambda, mu, xi_ani, phi_ani, eta_ani, &
if (anel_true) then
call plot_model_vtk(rho, lambda, mu, xi_ani, phi_ani, eta_ani, fa_ani_theta, &
- fa_ani_phi, Q_mu, Q_kappa)
+ fa_ani_phi, Q_mu_1d, Q_kappa_1d)
else
call plot_model_vtk(rho, lambda, mu, xi_ani, phi_ani, eta_ani, fa_ani_theta, &
fa_ani_phi)
@@ -272,6 +273,16 @@ subroutine read_model(rho, lambda, mu, xi_ani, phi_ani, eta_ani, &
if (use_netcdf) then
if (anel_true) then
+ do iel = 1, nelem ! iel
+ do ipol=0,npol ! ipol
+ do jpol=0,npol ! jpol
+ call compute_coordinates(s, z, r, theta, iel, ipol, jpol)
+ Q_mu(ipol,jpol,iel) = velocity(r, 'Qmu', iidom, modelstring, lfbkgrdmodel)
+ Q_kappa(ipol,jpol,iel) = velocity(r, 'Qka', iidom, modelstring, lfbkgrdmodel)
+ end do ! jpol
+ end do ! ipol
+ end do ! iel
+
call nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani, &
fa_ani_theta, fa_ani_phi, Q_mu, Q_kappa)
else
diff --git a/SOLVER/nc_routines.F90 b/SOLVER/nc_routines.F90
index 2bd0bd9..7837b27 100644
--- a/SOLVER/nc_routines.F90
+++ b/SOLVER/nc_routines.F90
@@ -52,6 +52,8 @@ module nc_routines
real(sp), allocatable :: scoord1d_mp(:), zcoord1d_mp(:)
real(sp), allocatable :: rho1d(:), mu1d(:), lambda1d(:)
real(sp), allocatable :: vp1d(:), vs1d(:)
+ real(sp), allocatable :: xi1d(:), phi1d(:), eta1d(:)
+ real(sp), allocatable :: Q_mu1d(:), Q_kappa1d(:)
!> Number of steps before kernel specific stuff is dumped
integer :: dumpstepsnap
@@ -766,7 +768,7 @@ subroutine nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani,
real(kind=dp), dimension(0:,0:,:), intent(in) :: rho, lambda, mu, xi_ani
real(kind=dp), dimension(0:,0:,:), intent(in) :: phi_ani, eta_ani
real(kind=dp), dimension(0:,0:,:), intent(in) :: fa_ani_theta, fa_ani_phi
- real(kind=sp), dimension(:), intent(in), optional :: Q_mu, Q_kappa
+ real(kind=dp), dimension(0:,0:,:), intent(in), optional :: Q_mu, Q_kappa
integer :: size1d
integer :: iel, ipol, jpol, ct
@@ -777,15 +779,29 @@ subroutine nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani,
allocate(mu1d(npoints))
allocate(vp1d(npoints))
allocate(vs1d(npoints))
+ allocate(xi1d(npoints))
+ allocate(phi1d(npoints))
+ allocate(eta1d(npoints))
+ if (present(Q_mu).and.present(Q_kappa)) then
+ allocate(Q_mu1d(npoints))
+ allocate(Q_kappa1d(npoints))
+ end if
do iel=1, nel_solid
do ipol=0, npol
do jpol=0, npol
if (kwf_mask(ipol,jpol,iel)) then
ct = mapping_ijel_ikwf(ipol,jpol,iel)
- rho1d(ct) = rho(ipol,jpol,ielsolid(iel))
+ rho1d(ct) = rho(ipol,jpol,ielsolid(iel))
lambda1d(ct) = lambda(ipol,jpol,ielsolid(iel))
- mu1d(ct) = mu(ipol,jpol,ielsolid(iel))
+ mu1d(ct) = mu(ipol,jpol,ielsolid(iel))
+ xi1d(ct) = xi_ani(ipol,jpol,ielsolid(iel))
+ phi1d(ct) = phi_ani(ipol,jpol,ielsolid(iel))
+ eta1d(ct) = eta_ani(ipol,jpol,ielsolid(iel))
+ if (present(Q_mu).and.present(Q_kappa)) then
+ Q_mu1d(ct) = Q_mu(ipol,jpol,ielsolid(iel))
+ Q_kappa1d(ct) = Q_kappa(ipol,jpol,ielsolid(iel))
+ endif
endif
enddo
enddo
@@ -796,9 +812,16 @@ subroutine nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani,
do jpol=0, npol
if (kwf_mask(ipol,jpol,iel + nel_solid)) then
ct = mapping_ijel_ikwf(ipol,jpol,iel + nel_solid)
- rho1d(ct) = rho(ipol,jpol,ielsolid(iel))
- lambda1d(ct) = lambda(ipol,jpol,ielsolid(iel))
- mu1d(ct) = mu(ipol,jpol,ielsolid(iel))
+ rho1d(ct) = rho(ipol,jpol,ielfluid(iel))
+ lambda1d(ct) = lambda(ipol,jpol,ielfluid(iel))
+ mu1d(ct) = mu(ipol,jpol,ielfluid(iel))
+ xi1d(ct) = xi_ani(ipol,jpol,ielfluid(iel))
+ phi1d(ct) = phi_ani(ipol,jpol,ielfluid(iel))
+ eta1d(ct) = eta_ani(ipol,jpol,ielfluid(iel))
+ if (present(Q_mu).and.present(Q_kappa)) then
+ Q_mu1d(ct) = Q_mu(ipol,jpol,ielfluid(iel))
+ Q_kappa1d(ct) = Q_kappa(ipol,jpol,ielfluid(iel))
+ endif
endif
enddo
enddo
@@ -815,12 +838,26 @@ subroutine nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani,
allocate(mu1d(size1d))
allocate(vp1d(size1d))
allocate(vs1d(size1d))
+ allocate(xi1d(size1d))
+ allocate(phi1d(size1d))
+ allocate(eta1d(size1d))
rho1d = real(pack(rho(ibeg:iend, ibeg:iend, :) ,.true.), kind=sp)
lambda1d = real(pack(lambda(ibeg:iend, ibeg:iend, :) ,.true.), kind=sp)
mu1d = real(pack(mu(ibeg:iend, ibeg:iend, :) ,.true.), kind=sp)
vp1d = sqrt( (lambda1d + 2.*mu1d ) / rho1d )
vs1d = sqrt( mu1d / rho1d )
+
+ xi1d = real(pack(xi_ani(ibeg:iend, ibeg:iend, :) ,.true.), kind=sp)
+ phi1d = real(pack(phi_ani(ibeg:iend, ibeg:iend, :) ,.true.), kind=sp)
+ eta1d = real(pack(eta_ani(ibeg:iend, ibeg:iend, :) ,.true.), kind=sp)
+
+ if (present(Q_mu).and.present(Q_kappa)) then
+ allocate(Q_mu1d(size1d))
+ allocate(Q_kappa1d(size1d))
+ Q_mu1d = real(pack(Q_mu(ibeg:iend, ibeg:iend, :) ,.true.), kind=sp)
+ Q_kappa1d = real(pack(Q_kappa(ibeg:iend, ibeg:iend, :) ,.true.), kind=sp)
+ endif
endif
end subroutine nc_dump_elastic_parameters
@@ -866,8 +903,7 @@ subroutine nc_define_outputfile(nrec, rec_names, rec_th, rec_th_req, rec_ph, rec
use data_io, only: datapath, lfdata, strain_samp
use data_mesh, only: maxind, num_rec, discont, nelem, nel_solid, nel_fluid, &
ndisc, maxind_glob, nelem_kwf_global, npoint_kwf, npoint_solid_kwf, &
- npoint_fluid_kwf, npol, nelem_kwf, npoint_kwf_global
-
+ npoint_fluid_kwf, npol, nelem_kwf, npoint_kwf_global, anel_true
use data_source, only: src_type, t_0
use data_time, only: deltat, niter
@@ -895,6 +931,11 @@ subroutine nc_define_outputfile(nrec, rec_names, rec_th, rec_th_req, rec_ph, rec
integer :: nc_mesh_vs_varid, nc_mesh_vp_varid
integer :: nc_mesh_mu_varid, nc_mesh_rho_varid
integer :: nc_mesh_lambda_varid
+ integer :: nc_mesh_xi_varid
+ integer :: nc_mesh_phi_varid
+ integer :: nc_mesh_eta_varid
+ integer :: nc_mesh_qmu_varid
+ integer :: nc_mesh_qka_varid
integer :: nc_mesh_midpoint_varid
integer :: nc_mesh_fem_varid
integer :: nc_mesh_eltype_varid
@@ -1227,6 +1268,35 @@ subroutine nc_define_outputfile(nrec, rec_names, rec_th, rec_th_req, rec_ph, rec
xtype = NF90_FLOAT, &
dimids = nc_pt_dimid,&
varid = nc_mesh_mu_varid) )
+ call check( nf90_def_var( ncid = ncid_meshout, &
+ name = 'mesh_xi', &
+ xtype = NF90_FLOAT, &
+ dimids = nc_pt_dimid,&
+ varid = nc_mesh_xi_varid) )
+ call check( nf90_def_var( ncid = ncid_meshout, &
+ name = 'mesh_phi', &
+ xtype = NF90_FLOAT, &
+ dimids = nc_pt_dimid,&
+ varid = nc_mesh_phi_varid) )
+ call check( nf90_def_var( ncid = ncid_meshout, &
+ name = 'mesh_eta', &
+ xtype = NF90_FLOAT, &
+ dimids = nc_pt_dimid,&
+ varid = nc_mesh_eta_varid) )
+
+ if (anel_true) then
+
+ call check( nf90_def_var( ncid = ncid_meshout, &
+ name = 'mesh_Qmu', &
+ xtype = NF90_FLOAT, &
+ dimids = nc_pt_dimid,&
+ varid = nc_mesh_Qmu_varid) )
+ call check( nf90_def_var( ncid = ncid_meshout, &
+ name = 'mesh_Qka', &
+ xtype = NF90_FLOAT, &
+ dimids = nc_pt_dimid,&
+ varid = nc_mesh_Qka_varid) )
+ end if
! call check( nf90_def_var( ncid = ncid_meshout, &
! name = 'model_domain', &
@@ -1524,6 +1594,11 @@ subroutine nc_finish_prepare
integer :: nc_mesh_vs_varid, nc_mesh_vp_varid
integer :: nc_mesh_mu_varid, nc_mesh_rho_varid
integer :: nc_mesh_lambda_varid
+ integer :: nc_mesh_xi_varid
+ integer :: nc_mesh_phi_varid
+ integer :: nc_mesh_eta_varid
+ integer :: nc_mesh_qmu_varid
+ integer :: nc_mesh_qka_varid
integer :: nc_mesh_midpoint_varid
integer :: nc_mesh_eltype_varid
@@ -1648,6 +1723,50 @@ subroutine nc_finish_prepare
start = npoints_myfirst, &
count = npoints )
+ ! Anisotropic parameters
+ ! Phi
+ call getvarid( ncid_meshout, "mesh_phi", nc_mesh_phi_varid )
+ call putvar_real1d( ncid = ncid_meshout, &
+ varid = nc_mesh_mu_varid, &
+ values = phi1d, &
+ start = npoints_myfirst, &
+ count = npoints )
+
+ ! Xi
+ call getvarid( ncid_meshout, "mesh_xi", nc_mesh_xi_varid )
+ call putvar_real1d( ncid = ncid_meshout, &
+ varid = nc_mesh_xi_varid, &
+ values = xi1d, &
+ start = npoints_myfirst, &
+ count = npoints )
+
+ ! Eta
+ call getvarid( ncid_meshout, "mesh_eta", nc_mesh_eta_varid )
+ call putvar_real1d( ncid = ncid_meshout, &
+ varid = nc_mesh_eta_varid, &
+ values = eta1d, &
+ start = npoints_myfirst, &
+ count = npoints )
+
+ ! Anelastic parameters
+ if (allocated(Q_mu1d).and.allocated(Q_kappa1d)) then
+
+ ! Q_mu
+ call getvarid( ncid_meshout, "mesh_Qmu", nc_mesh_Qmu_varid )
+ call putvar_real1d( ncid = ncid_meshout, &
+ varid = nc_mesh_Qmu_varid, &
+ values = Q_mu1d, &
+ start = npoints_myfirst, &
+ count = npoints )
+ ! Q_kappa
+ call getvarid( ncid_meshout, "mesh_Qka", nc_mesh_Qka_varid )
+ call putvar_real1d( ncid = ncid_meshout, &
+ varid = nc_mesh_Qka_varid, &
+ values = Q_kappa1d, &
+ start = npoints_myfirst, &
+ count = npoints )
+ end if
+
if (trim(dump_type) == 'displ_only' .and. nelem_kwf > 0) then
call getvarid( ncid_meshout, "midpoint_mesh", nc_mesh_midpoint_varid )
call check(nf90_put_var ( ncid = ncid_meshout, &
More information about the CIG-COMMITS
mailing list