[cig-commits] [commit] master: (re)introducing jbeg and jend for kernel wavefields in (7a0ed41)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Oct 13 09:07:42 PDT 2014
Repository : https://github.com/geodynamics/axisem
On branch : master
Link : https://github.com/geodynamics/axisem/compare/ca4b3002407f17f95bbf7013f7622a43ceb4db35...7a0ed4115e0445ccb0219c894eedf4bf157d8689
>---------------------------------------------------------------
commit 7a0ed4115e0445ccb0219c894eedf4bf157d8689
Author: martinvandriel <martin at vandriel.de>
Date: Mon Oct 13 18:00:57 2014 +0200
(re)introducing jbeg and jend for kernel wavefields in
fullfields or strain_only mode. Useful for instaseis to get denser lateral then depth
sampling
>---------------------------------------------------------------
7a0ed4115e0445ccb0219c894eedf4bf157d8689
SOLVER/data_io.f90 | 5 ++--
SOLVER/get_mesh.f90 | 4 ++-
SOLVER/get_model.F90 | 2 +-
SOLVER/inparam_advanced.TEMPLATE | 2 ++
SOLVER/meshes_io.F90 | 45 ++++++++++++++--------------
SOLVER/nc_routines.F90 | 26 ++++++++--------
SOLVER/parameters.F90 | 20 +++++++++++--
SOLVER/time_evol_wave.F90 | 4 +--
SOLVER/wavefields_io.f90 | 64 ++++++++++++++++++++--------------------
9 files changed, 96 insertions(+), 76 deletions(-)
diff --git a/SOLVER/data_io.f90 b/SOLVER/data_io.f90
index 6eea9b3..038a151 100644
--- a/SOLVER/data_io.f90
+++ b/SOLVER/data_io.f90
@@ -61,8 +61,9 @@ module data_io
! indices to limit dumping to select contiguous range of GLL points:
! 0<=ibeg<=iend<=npol
! For the time being: dump the same in xeta and eta directions
- integer :: ibeg,iend
- ! ndumppts_el=(iend-ibeg+1)**2
+ integer :: ibeg, iend
+ integer :: jbeg, jend
+ ! ndumppts_el = (iend-ibeg+1) * (jend-jbeg+1)
integer :: ndumppts_el
! for xdmf dumps
diff --git a/SOLVER/get_mesh.f90 b/SOLVER/get_mesh.f90
index 32e5cb1..fb25352 100644
--- a/SOLVER/get_mesh.f90
+++ b/SOLVER/get_mesh.f90
@@ -54,7 +54,7 @@ subroutine read_db
use data_comm
use data_proc
use data_time
- use data_io, only : do_anel, ibeg, iend, dump_type
+ use data_io, only : do_anel, ibeg, iend, jbeg, jend, dump_type
use commun, only : barrier, psum, pmax, pmin
use background_models, only : model_is_ani, model_is_anelastic, get_ext_disc, &
override_ext_q
@@ -92,6 +92,8 @@ subroutine read_db
if (trim(dump_type) == 'displ_only') then
ibeg = 0
iend = npol
+ jbeg = 0
+ jend = npol
endif
!!!!!!!!!!!! BACKGROUND MODEL !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
diff --git a/SOLVER/get_model.F90 b/SOLVER/get_model.F90
index 692baaf..aa72b8b 100644
--- a/SOLVER/get_model.F90
+++ b/SOLVER/get_model.F90
@@ -226,7 +226,7 @@ subroutine read_model(rho, lambda, mu, xi_ani, phi_ani, eta_ani, &
! write out for later snaps
if (do_mesh_tests) then
do ipol=ibeg, iend
- do jpol=ibeg, iend
+ do jpol=jbeg, jend
call compute_coordinates(s,z,r,theta,iel,ipol,jpol)
write(60000+mynum,14)r,theta,rho(ipol,jpol,iel),&
sqrt( (lambda(ipol,jpol,iel)+2.*mu(ipol,jpol,iel))/rho(ipol,jpol,iel)), &
diff --git a/SOLVER/inparam_advanced.TEMPLATE b/SOLVER/inparam_advanced.TEMPLATE
index 15fbb98..02ee7c5 100644
--- a/SOLVER/inparam_advanced.TEMPLATE
+++ b/SOLVER/inparam_advanced.TEMPLATE
@@ -141,6 +141,8 @@ KERNEL_SOURCE mask
# (overwritten with 0 and npol for dumptype displ_only)
KERNEL_IBEG 1
KERNEL_IEND 3
+KERNEL_JBEG 1
+KERNEL_JEND 3
# minimal and maximal colatitude in degree for kernel wavefields
# (only for dumptype displ_only)
diff --git a/SOLVER/meshes_io.F90 b/SOLVER/meshes_io.F90
index 69c7917..5ef64e3 100644
--- a/SOLVER/meshes_io.F90
+++ b/SOLVER/meshes_io.F90
@@ -492,7 +492,7 @@ subroutine build_kwf_grid()
use nc_routines, only : set_npoints
use data_mesh
- use data_io, only: ibeg, iend
+ use data_io, only: ibeg, iend, jbeg, jend
integer :: iel, ipol, jpol, ct, ipt, idest
integer, allocatable :: mapping(:)
@@ -556,7 +556,7 @@ subroutine build_kwf_grid()
do iel=1, nel_solid
if (.not. mask_tp_elem(iel)) cycle
- do jpol=ibeg, iend
+ do jpol=jbeg, jend
do ipol=ibeg, iend
ipt = (iel-1)*(npol+1)**2 + jpol*(npol+1) + ipol + 1
@@ -595,7 +595,7 @@ subroutine build_kwf_grid()
do iel=1, nel_fluid
if (.not. mask_tp_elem(iel + nel_solid)) cycle
- do jpol=ibeg, iend
+ do jpol=jbeg, jend
do ipol=ibeg, iend
ipt = (iel-1)*(npol+1)**2 + jpol*(npol+1) + ipol + 1
@@ -635,7 +635,7 @@ subroutine build_kwf_grid()
call set_npoints(npoint_kwf)
if (lpr) then
- write(6,*) 'local point number: ', nelem_kwf * (iend - ibeg + 1)**2
+ write(6,*) 'local point number: ', nelem_kwf * (iend - ibeg + 1) * (jend - jbeg + 1)
write(6,*) 'after removing duplicates: ', npoint_kwf
write(6,*) 'compression: ', &
real(npoint_kwf) / real(nelem_kwf * (npol + 1)**2)
@@ -749,7 +749,7 @@ subroutine build_kwf_grid()
enddo
endif
- allocate(sem_mesh_kwf(ibeg:iend, ibeg:iend, 1:nelem_kwf))
+ allocate(sem_mesh_kwf(ibeg:iend, jbeg:jend, 1:nelem_kwf))
if (lpr) write(6,*) ' .... constructing spectral element grid for kwf output'
@@ -758,7 +758,7 @@ subroutine build_kwf_grid()
do iel=1, nel_solid
if (.not. mask_tp_elem(iel)) cycle
do ipol=ibeg, iend
- do jpol=ibeg, iend
+ do jpol=jbeg, jend
sem_mesh_kwf(ipol,jpol,ct) = mapping_ijel_ikwf(ipol,jpol,iel) - 1
enddo
enddo
@@ -768,7 +768,7 @@ subroutine build_kwf_grid()
do iel=1, nel_fluid
if (.not. mask_tp_elem(iel + nel_solid)) cycle
do ipol=ibeg, iend
- do jpol=ibeg, iend
+ do jpol=jbeg, jend
sem_mesh_kwf(ipol,jpol,ct) = mapping_ijel_ikwf(ipol,jpol,iel + nel_solid) - 1
enddo
enddo
@@ -1016,7 +1016,7 @@ subroutine dump_kwf_sem_xdmf(filename, npoints, nelem)
open(newunit=iinput_xdmf, file=trim(filename)//'_sem.xdmf')
write(iinput_xdmf, 733) npoints, npoints, trim(filename_np), npoints, trim(filename_np)
- write(iinput_xdmf, 734) nelem, (iend-ibeg+1)**2, nelem, iend-ibeg+1, iend-ibeg+1, &
+ write(iinput_xdmf, 734) nelem, (npol+1)**2, nelem, npol+1, npol+1, &
trim(filename_np), "'", "'"
close(iinput_xdmf)
@@ -1188,7 +1188,7 @@ end subroutine
subroutine dump_wavefields_mesh_1d
use data_mesh
- use data_io, only: ibeg, iend, ndumppts_el
+ use data_io, only: ibeg, iend, jbeg, jend, ndumppts_el
use data_spec, only: G1T, G2T, G2
use data_pointwise
use nc_routines, only: nc_dump_mesh_sol, nc_dump_mesh_flu
@@ -1203,13 +1203,14 @@ subroutine dump_wavefields_mesh_1d
if (lpr) then
write(6,*)' set strain dumping GLL boundaries to:'
- write(6,*)' ipol=',ibeg,iend
+ write(6,*)' ipol=', ibeg, iend
+ write(6,*)' jpol=', jbeg, jend
endif
- allocate(ssol(ibeg:iend,ibeg:iend,nel_solid))
- allocate(zsol(ibeg:iend,ibeg:iend,nel_solid))
- allocate(sflu(ibeg:iend,ibeg:iend,nel_fluid))
- allocate(zflu(ibeg:iend,ibeg:iend,nel_fluid))
+ allocate(ssol(ibeg:iend,jbeg:jend,nel_solid))
+ allocate(zsol(ibeg:iend,jbeg:jend,nel_solid))
+ allocate(sflu(ibeg:iend,jbeg:jend,nel_fluid))
+ allocate(zflu(ibeg:iend,jbeg:jend,nel_fluid))
endif
@@ -1218,7 +1219,7 @@ subroutine dump_wavefields_mesh_1d
else
! compute solid grid
do iel=1,nel_solid
- do jpol=ibeg,iend
+ do jpol=jbeg,jend
do ipol=ibeg,iend
ssol(ipol,jpol,iel) = scoord(ipol,jpol,ielsolid(iel))
zsol(ipol,jpol,iel) = zcoord(ipol,jpol,ielsolid(iel))
@@ -1228,15 +1229,15 @@ subroutine dump_wavefields_mesh_1d
if (lpr) write(6,*)' dumping solid submesh for kernel wavefields...'
if (use_netcdf) then
- call nc_dump_mesh_sol(real(ssol(ibeg:iend,ibeg:iend,:)), &
- real(zsol(ibeg:iend,ibeg:iend,:)))
+ call nc_dump_mesh_sol(real(ssol(ibeg:iend,jbeg:jend,:)), &
+ real(zsol(ibeg:iend,jbeg:jend,:)))
else
open(unit=2500+mynum,file=datapath(1:lfdata)//'/strain_mesh_sol_'&
//appmynum//'.dat', &
FORM="UNFORMATTED",STATUS="REPLACE")
- write(2500+mynum)ssol(ibeg:iend,ibeg:iend,:),zsol(ibeg:iend,ibeg:iend,:)
+ write(2500+mynum)ssol(ibeg:iend,jbeg:jend,:),zsol(ibeg:iend,jbeg:jend,:)
close(2500+mynum)
end if
deallocate(ssol,zsol)
@@ -1244,7 +1245,7 @@ subroutine dump_wavefields_mesh_1d
! compute fluid grid
if (have_fluid) then
do iel=1,nel_fluid
- do jpol=ibeg,iend
+ do jpol=jbeg,jend
do ipol=ibeg,iend
sflu(ipol,jpol,iel) = scoord(ipol,jpol,ielfluid(iel))
zflu(ipol,jpol,iel) = zcoord(ipol,jpol,ielfluid(iel))
@@ -1253,13 +1254,13 @@ subroutine dump_wavefields_mesh_1d
enddo
if (lpr) write(6,*)' dumping fluid submesh for kernel wavefields...'
if (use_netcdf) then
- call nc_dump_mesh_flu(real(sflu(ibeg:iend,ibeg:iend,:)),&
- real(zflu(ibeg:iend,ibeg:iend,:)))
+ call nc_dump_mesh_flu(real(sflu(ibeg:iend,jbeg:jend,:)),&
+ real(zflu(ibeg:iend,jbeg:jend,:)))
else
open(unit=2600+mynum,file=datapath(1:lfdata)//'/strain_mesh_flu_'&
//appmynum//'.dat', &
FORM="UNFORMATTED",STATUS="REPLACE")
- write(2600+mynum)sflu(ibeg:iend,ibeg:iend,:),zflu(ibeg:iend,ibeg:iend,:)
+ write(2600+mynum)sflu(ibeg:iend,jbeg:jend,:),zflu(ibeg:iend,jbeg:jend,:)
close(2600+mynum)
end if
deallocate(sflu,zflu)
diff --git a/SOLVER/nc_routines.F90 b/SOLVER/nc_routines.F90
index 6f46ee9..7a42576 100644
--- a/SOLVER/nc_routines.F90
+++ b/SOLVER/nc_routines.F90
@@ -765,7 +765,7 @@ end subroutine nc_dump_mesh_mp_kwf
subroutine nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani, &
fa_ani_theta, fa_ani_phi, Q_mu, Q_kappa)
- use data_io, only: ibeg, iend, dump_type
+ use data_io, only: ibeg, iend, jbeg, jend, dump_type
use data_mesh, only: mapping_ijel_ikwf, ielsolid, ielfluid, nel_solid, nel_fluid, &
npol, kwf_mask
@@ -835,7 +835,7 @@ subroutine nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani,
vs1d = sqrt( mu1d / rho1d )
else
- size1d = size(rho(ibeg:iend, ibeg:iend, :))
+ size1d = size(rho(ibeg:iend, jbeg:jend, :))
print *, ' NetCDF: Mesh elastic parameter variables have size:', size1d
allocate(rho1d(size1d))
allocate(lambda1d(size1d))
@@ -846,21 +846,21 @@ subroutine nc_dump_elastic_parameters(rho, lambda, mu, xi_ani, phi_ani, eta_ani,
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)
+ rho1d = real(pack(rho(ibeg:iend, jbeg:jend, :) ,.true.), kind=sp)
+ lambda1d = real(pack(lambda(ibeg:iend, jbeg:jend, :) ,.true.), kind=sp)
+ mu1d = real(pack(mu(ibeg:iend, jbeg:jend, :) ,.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)
+ xi1d = real(pack(xi_ani(ibeg:iend, jbeg:jend, :) ,.true.), kind=sp)
+ phi1d = real(pack(phi_ani(ibeg:iend, jbeg:jend, :) ,.true.), kind=sp)
+ eta1d = real(pack(eta_ani(ibeg:iend, jbeg:jend, :) ,.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)
+ Q_mu1d = real(pack(Q_mu(ibeg:iend, jbeg:jend, :) ,.true.), kind=sp)
+ Q_kappa1d = real(pack(Q_kappa(ibeg:iend, jbeg:jend, :) ,.true.), kind=sp)
endif
endif
@@ -902,8 +902,8 @@ end subroutine nc_dump_elastic_parameters
!! and allocate buffer variables.
subroutine nc_define_outputfile(nrec, rec_names, rec_th, rec_th_req, rec_ph, rec_proc)
- use data_io, only: nseismo, nstrain, nseismo, ibeg, iend, dump_wavefields, &
- dump_type
+ use data_io, only: nseismo, nstrain, nseismo, ibeg, iend, jbeg, jend, &
+ dump_wavefields, dump_type
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, &
@@ -1153,7 +1153,7 @@ subroutine nc_define_outputfile(nrec, rec_names, rec_th, rec_th_req, rec_ph, rec
'velo_s ', 'velo_p ', 'velo_z ']
end if
- gllperelem = (iend - ibeg + 1)**2
+ gllperelem = (iend - ibeg + 1) * (jend - jbeg + 1)
npoints = nelem * gllperelem
call comm_elem_number(npoints, npoints_global, npoints_myfirst, npoints_mylast)
diff --git a/SOLVER/parameters.F90 b/SOLVER/parameters.F90
index 5552beb..e65ae40 100644
--- a/SOLVER/parameters.F90
+++ b/SOLVER/parameters.F90
@@ -99,7 +99,7 @@ subroutine readin_parameters
enforced_period, trim(simtype), rec_file_type, &
sum_seis, sum_fields, time_scheme, seis_dt, &
dump_energy, dump_vtk, dump_wavefields, &
- dump_type, ibeg, iend, strain_samp, src_dump_type, make_homo, &
+ dump_type, ibeg, iend, jbeg, jend, strain_samp, src_dump_type, make_homo, &
add_hetero, do_mesh_tests, output_format
20 format(/&
@@ -198,6 +198,8 @@ subroutine readin_parameters
12x,'Wavefield dumping type: ',a12,/ &
12x,'First GLL to save in strains: ',i2,/ &
12x,'Last GLL to save in strains: ',i2,/ &
+ 12x,'First GLL to save in strains: ',i2,/ &
+ 12x,'Last GLL to save in strains: ',i2,/ &
12x,'Samples per period for strains: ',f7.3,/ &
12x,'Source dumping type: ',a4,/ &
12x,'Homogenize background model? ',l2,/ &
@@ -384,6 +386,8 @@ subroutine read_inparam_advanced
src_dump_type = 'mask'
ibeg = 1
iend = 3
+ jbeg = 1
+ jend = 3
kwf_rmin = 0
kwf_rmax = 7d6
@@ -488,7 +492,12 @@ subroutine read_inparam_advanced
case('KERNEL_IEND')
read(keyvalue,*) iend
- !iend = npol - iend
+
+ case('KERNEL_JBEG')
+ read(keyvalue,*) jbeg
+
+ case('KERNEL_JEND')
+ read(keyvalue,*) jend
case('KERNEL_RMIN')
read(keyvalue, *) kwf_rmin
@@ -609,6 +618,8 @@ subroutine read_inparam_advanced
call broadcast_int(ibeg, 0)
call broadcast_int(iend, 0)
+ call broadcast_int(jbeg, 0)
+ call broadcast_int(jend, 0)
call broadcast_log(dump_energy, 0)
call broadcast_log(make_homo, 0)
@@ -1088,10 +1099,11 @@ subroutine compute_numerical_parameters
write(6,13)' ...that is, every :',strain_it,'timestep'
endif
- ndumppts_el = (iend - ibeg + 1)**2
+ ndumppts_el = (iend - ibeg + 1) * (jend - jbeg + 1)
if (lpr) then
write(6,*)' Define limitation of GLL points in the dumped fields:'
write(6,*)' ibeg=', ibeg, 'iend=', iend
+ write(6,*)' jbeg=', jbeg, 'jend=', jend
write(6,*)' # points saved within an element:', ndumppts_el
endif
endif
@@ -1501,6 +1513,8 @@ subroutine write_parameters
call nc_write_att_char( 'cyl', 'receiver components ')
call nc_write_att_int( ibeg, 'ibeg')
call nc_write_att_int( iend, 'iend')
+ call nc_write_att_int( jbeg, 'jbeg')
+ call nc_write_att_int( jend, 'jend')
call nc_write_att_real( shift_fact, 'source shift factor in sec')
call nc_write_att_int( nint(shift_fact/deltat), 'source shift factor for deltat')
call nc_write_att_int( nint(shift_fact/seis_dt), 'source shift factor for seis_dt')
diff --git a/SOLVER/time_evol_wave.F90 b/SOLVER/time_evol_wave.F90
index 55d5acb..1c96d08 100644
--- a/SOLVER/time_evol_wave.F90
+++ b/SOLVER/time_evol_wave.F90
@@ -111,7 +111,7 @@ subroutine prepare_waves
if (dump_vtk) then
if (lpr) write(6,*)' dumping global grids for snapshots...'
- call dump_glob_grid_midpoint(ibeg,iend,ibeg,iend)
+ call dump_glob_grid_midpoint(ibeg,iend,jbeg,jend)
endif
if (dump_xdmf) then
@@ -1081,7 +1081,7 @@ subroutine dump_stuff(iter, iseismo, istrain, isnap, &
write(6,*) 'Writing global snap to file: ', isnap
write(6,*)
endif
- call glob_snapshot_midpoint(disp, chi, ibeg, iend, ibeg, iend, isnap)
+ call glob_snapshot_midpoint(disp, chi, ibeg, iend, jbeg, jend, isnap)
endif
endif
diff --git a/SOLVER/wavefields_io.f90 b/SOLVER/wavefields_io.f90
index 3341fb1..6ed92ab 100644
--- a/SOLVER/wavefields_io.f90
+++ b/SOLVER/wavefields_io.f90
@@ -790,7 +790,7 @@ subroutine dump_field_1d(f, filename, appisnap, n)
if (npoint_solid_kwf > 0) &
call nc_dump_field_solid(kwf_mapping_sol(floc), filename(2:))
else
- call nc_dump_field_solid(pack(floc(ibeg:iend,ibeg:iend,:), .true.), &
+ call nc_dump_field_solid(pack(floc(ibeg:iend,jbeg:jend,:), .true.), &
filename(2:))
endif
@@ -799,7 +799,7 @@ subroutine dump_field_1d(f, filename, appisnap, n)
if (npoint_fluid_kwf > 0) &
call nc_dump_field_fluid(kwf_mapping_flu(floc), filename(2:))
else
- call nc_dump_field_fluid(pack(floc(ibeg:iend,ibeg:iend,:), .true.), &
+ call nc_dump_field_fluid(pack(floc(ibeg:iend,jbeg:jend,:), .true.), &
filename(2:))
endif
@@ -811,7 +811,7 @@ subroutine dump_field_1d(f, filename, appisnap, n)
open(unit=25000+mynum, file=datapath(1:lfdata)//filename//'_' &
//appmynum//'_'//appisnap//'.bindat', &
FORM="UNFORMATTED", STATUS="UNKNOWN", POSITION="REWIND")
- write(25000+mynum) pack(floc(ibeg:iend,ibeg:iend,:), .true.)
+ write(25000+mynum) pack(floc(ibeg:iend,jbeg:jend,:), .true.)
close(25000+mynum)
end if
@@ -845,9 +845,9 @@ subroutine dump_disp(u, chi, istrain)
FORM="UNFORMATTED",STATUS="REPLACE")
if (src_type(1)/='monopole') then
- write(75000+mynum) (f(ibeg:iend,ibeg:iend,:,i),i=1,3)
+ write(75000+mynum) (f(ibeg:iend,jbeg:jend,:,i),i=1,3)
else
- write(75000+mynum) f(ibeg:iend,ibeg:iend,:,1:3:2)
+ write(75000+mynum) f(ibeg:iend,jbeg:jend,:,1:3:2)
endif
close(75000+mynum)
@@ -891,9 +891,9 @@ subroutine dump_velo_dchi(v, dchi, istrain)
FORM="UNFORMATTED",STATUS="REPLACE")
if (src_type(1)/='monopole') then
- write(85000+mynum) (f(ibeg:iend,ibeg:iend,:,i), i=1,3)
+ write(85000+mynum) (f(ibeg:iend,jbeg:jend,:,i), i=1,3)
else
- write(85000+mynum) f(ibeg:iend,ibeg:iend,:,1), f(ibeg:iend,ibeg:iend,:,3)
+ write(85000+mynum) f(ibeg:iend,jbeg:jend,:,1), f(ibeg:iend,jbeg:jend,:,3)
endif
close(85000+mynum)
@@ -947,19 +947,19 @@ subroutine dump_velo_global(v, dchi, istrain)
if (use_netcdf) then
if (src_type(1)/='monopole') then
- call nc_dump_field_solid(pack(f(ibeg:iend,ibeg:iend,:,2),.true.), 'velo_sol_p')
+ call nc_dump_field_solid(pack(f(ibeg:iend,jbeg:jend,:,2),.true.), 'velo_sol_p')
end if
- call nc_dump_field_solid(pack(f(ibeg:iend,ibeg:iend,:,1),.true.), 'velo_sol_s')
- call nc_dump_field_solid(pack(f(ibeg:iend,ibeg:iend,:,3),.true.), 'velo_sol_z')
+ call nc_dump_field_solid(pack(f(ibeg:iend,jbeg:jend,:,1),.true.), 'velo_sol_s')
+ call nc_dump_field_solid(pack(f(ibeg:iend,jbeg:jend,:,3),.true.), 'velo_sol_z')
else
open(unit=95000+mynum,file=datapath(1:lfdata)//'/velo_sol_'&
//appmynum//'_'//appisnap//'.bindat',&
FORM="UNFORMATTED",STATUS="REPLACE")
if (src_type(1)/='monopole') then
- write(95000+mynum) (f(ibeg:iend,ibeg:iend,:,i), i=1,3)
+ write(95000+mynum) (f(ibeg:iend,jbeg:jend,:,i), i=1,3)
else
- write(95000+mynum) f(ibeg:iend,ibeg:iend,:,1), &
- f(ibeg:iend,ibeg:iend,:,3)
+ write(95000+mynum) f(ibeg:iend,jbeg:jend,:,1), &
+ f(ibeg:iend,jbeg:jend,:,3)
end if
close(95000+mynum)
end if
@@ -971,25 +971,25 @@ subroutine dump_velo_global(v, dchi, istrain)
call axisym_gradient_fluid(dchi, usz_fluid)
call define_io_appendix(appisnap,istrain)
- fflu(ibeg:iend,ibeg:iend,:,1) = inv_rho_fluid(ibeg:iend,ibeg:iend,:) * &
- usz_fluid(ibeg:iend,ibeg:iend,:,1)
- fflu(ibeg:iend,ibeg:iend,:,2) = 0
- fflu(ibeg:iend,ibeg:iend,:,3) = inv_rho_fluid(ibeg:iend,ibeg:iend,:) * &
- usz_fluid(ibeg:iend,ibeg:iend,:,2)
+ fflu(ibeg:iend,jbeg:jend,:,1) = inv_rho_fluid(ibeg:iend,jbeg:jend,:) * &
+ usz_fluid(jbeg:jend,jbeg:jend,:,1)
+ fflu(ibeg:iend,jbeg:jend,:,2) = 0
+ fflu(ibeg:iend,jbeg:jend,:,3) = inv_rho_fluid(ibeg:iend,jbeg:jend,:) * &
+ usz_fluid(ibeg:iend,jbeg:jend,:,2)
! dump velocity vector inside fluid
if (use_netcdf) then
- call nc_dump_field_fluid(pack(fflu(ibeg:iend,ibeg:iend,:,1), .true.), 'velo_flu_s')
- call nc_dump_field_fluid(pack(fflu(ibeg:iend,ibeg:iend,:,3), .true.), 'velo_flu_z')
+ call nc_dump_field_fluid(pack(fflu(ibeg:iend,jbeg:jend,:,1), .true.), 'velo_flu_s')
+ call nc_dump_field_fluid(pack(fflu(ibeg:iend,jbeg:jend,:,3), .true.), 'velo_flu_z')
if (src_type(1)/='monopole') then
- call nc_dump_field_fluid(pack(fflu(ibeg:iend,ibeg:iend,:,2), .true.), 'velo_flu_p')
+ call nc_dump_field_fluid(pack(fflu(ibeg:iend,jbeg:jend,:,2), .true.), 'velo_flu_p')
end if
else
open(unit=960000+mynum,file=datapath(1:lfdata)//'/velo_flu_'&
//appmynum//'_'//appisnap//'.bindat',&
FORM="UNFORMATTED",STATUS="REPLACE")
- write(960000+mynum) (fflu(ibeg:iend,ibeg:iend,:,i), i=1,3)
+ write(960000+mynum) (fflu(ibeg:iend,jbeg:jend,:,i), i=1,3)
close(960000+mynum)
end if ! netcdf
endif ! have_fluid
@@ -1041,20 +1041,20 @@ subroutine dump_disp_global(u, chi, istrain)
elseif (use_netcdf .and. dump_type /= 'displ_only') then
if (src_type(1)/='monopole') then
- call nc_dump_field_solid(pack(f(ibeg:iend,ibeg:iend,:,2),.true.), 'disp_sol_p')
+ call nc_dump_field_solid(pack(f(ibeg:iend,jbeg:jend,:,2),.true.), 'disp_sol_p')
end if
- call nc_dump_field_solid(pack(f(ibeg:iend,ibeg:iend,:,1),.true.), 'disp_sol_s')
- call nc_dump_field_solid(pack(f(ibeg:iend,ibeg:iend,:,3),.true.), 'disp_sol_z')
+ call nc_dump_field_solid(pack(f(ibeg:iend,jbeg:jend,:,1),.true.), 'disp_sol_s')
+ call nc_dump_field_solid(pack(f(ibeg:iend,jbeg:jend,:,3),.true.), 'disp_sol_z')
else
open(unit=95000+mynum,file=datapath(1:lfdata)//'/disp_sol_'&
//appmynum//'_'//appisnap//'.bindat',&
FORM="UNFORMATTED",STATUS="REPLACE")
if (src_type(1)/='monopole') then
- write(95000+mynum) (f(ibeg:iend,ibeg:iend,:,i), i=1,3)
+ write(95000+mynum) (f(ibeg:iend,jbeg:jend,:,i), i=1,3)
else
- write(95000+mynum) f(ibeg:iend,ibeg:iend,:,1), &
- f(ibeg:iend,ibeg:iend,:,3)
+ write(95000+mynum) f(ibeg:iend,jbeg:jend,:,1), &
+ f(ibeg:iend,jbeg:jend,:,3)
end if
close(95000+mynum)
end if
@@ -1079,17 +1079,17 @@ subroutine dump_disp_global(u, chi, istrain)
call nc_dump_field_fluid(kwf_mapping_flu(fflu(:,:,:,3)), 'disp_flu_z')
elseif (use_netcdf .and. dump_type /= 'displ_only') then
- call nc_dump_field_fluid(pack(fflu(ibeg:iend,ibeg:iend,:,1), .true.), 'disp_flu_s')
- call nc_dump_field_fluid(pack(fflu(ibeg:iend,ibeg:iend,:,3), .true.), 'disp_flu_z')
+ call nc_dump_field_fluid(pack(fflu(ibeg:iend,jbeg:jend,:,1), .true.), 'disp_flu_s')
+ call nc_dump_field_fluid(pack(fflu(ibeg:iend,jbeg:jend,:,3), .true.), 'disp_flu_z')
if (src_type(1)/='monopole') then
- call nc_dump_field_fluid(pack(fflu(ibeg:iend,ibeg:iend,:,2), .true.), 'disp_flu_p')
+ call nc_dump_field_fluid(pack(fflu(ibeg:iend,jbeg:jend,:,2), .true.), 'disp_flu_p')
end if
else
open(unit=960000+mynum,file=datapath(1:lfdata)//'/disp_flu_'&
//appmynum//'_'//appisnap//'.bindat',&
FORM="UNFORMATTED",STATUS="REPLACE")
- write(960000+mynum) (fflu(ibeg:iend,ibeg:iend,:,i), i=1,3)
+ write(960000+mynum) (fflu(ibeg:iend,jbeg:jend,:,i), i=1,3)
close(960000+mynum)
end if ! netcdf
endif ! have_fluid
More information about the CIG-COMMITS
mailing list