Merge branch 'master' into release_1.1

Pull intermediate updates of master branch into release branch. In case that should ever be released, it might be useful to have the latest version.

Conflicts:
	MANUAL/manual_axisem1.1.tex

Author: Simon Stähler

SOLVER/UTILS/field_transform.F90

@@ -44,7 +44,7 @@ program field_transformation
     integer                         :: ncout_mesh_mpz_varid, ncin_mesh_mpz_varid
     integer                         :: ncout_surf_grpid, ncout_surfstrain_dimid, ncout_surf_dimid
     integer                         :: ncout_comp_dimid
-    integer                         :: ncout_surf_varids(6), ncin_surf_grpid, ncin_surf_varids(6)
+    integer                         :: ncout_surf_varids(7), ncin_surf_grpid, ncin_surf_varids(7)
 
     integer                         :: nc_mesh_npol_dimid, nc_mesh_cntrlpts_dimid, &
                                        nc_mesh_elem_dimid
@@ -341,9 +341,10 @@ program field_transformation
                               len    = nsurfelem,        &
                               dimid  = ncout_surf_dimid) )     
 
-    allocate(varname_surf(6))
+    allocate(varname_surf(7))
     varname_surf = ['elem_theta      ', 'displacement    ', 'velocity        ', &
-                    'disp_src        ', 'strain          ', 'stf_dump        ']
+                    'disp_src        ', 'strain          ', 'stf_dump        ', &
+                    'stf_d_dump      ']
     do ivar = 1, size(varname_surf)
         call check( nf90_inq_varid(ncid  = ncin_surf_grpid,        &
                                    varid = ncin_surf_varids(ivar), & 
@@ -385,16 +386,18 @@ program field_transformation
                                       shuffle = 1, deflate = 1,          &
                                       deflate_level = deflate_level) )
 
-    call check( nf90_def_var(ncid       = ncout_surf_grpid,   &
-                             name       = varname_surf(6),    &
-                             xtype      = NF90_FLOAT,         &
-                             dimids     = [ncout_snap_dimid], &
-                             chunksizes = [nsnap],            &
-                             varid      = ncout_surf_varids(6)) )
-    call check( nf90_def_var_deflate( ncid    = ncout_surf_grpid,        &
-                                      varid   = ncout_surf_varids(6),    &
-                                      shuffle = 1, deflate = 1,          &
-                                      deflate_level = deflate_level) )
+    do ivar = 6, 7
+       call check( nf90_def_var(ncid       = ncout_surf_grpid,   &
+                                name       = varname_surf(ivar),    &
+                                xtype      = NF90_FLOAT,         &
+                                dimids     = [ncout_snap_dimid], &
+                                chunksizes = [nsnap],            &
+                                varid      = ncout_surf_varids(ivar)) )
+       call check( nf90_def_var_deflate( ncid    = ncout_surf_grpid,        &
+                                         varid   = ncout_surf_varids(ivar),    &
+                                         shuffle = 1, deflate = 1,          &
+                                         deflate_level = deflate_level) )
+    enddo
 
     print *, 'Surface variables defined'
 
@@ -608,16 +611,18 @@ program field_transformation
     deallocate(data_surf_3d)
 
     allocate(data_surf_1d(nsnap))
-    call check( nf90_get_var( ncid   = ncin_surf_grpid,        &
-                              varid  = ncin_surf_varids(6),    &
-                              start  = [1],                    & 
-                              count  = [nsnap],                &
-                              values = data_surf_1d) )
-    call check( nf90_put_var( ncid   = ncout_surf_grpid,       &
-                              varid  = ncout_surf_varids(6),   &
-                              start  = [1],                    & 
-                              count  = [nsnap],                &
-                              values = data_surf_1d))
+    do ivar = 6, 7
+       call check( nf90_get_var( ncid   = ncin_surf_grpid,        &
+                                 varid  = ncin_surf_varids(ivar),    &
+                                 start  = [1],                    & 
+                                 count  = [nsnap],                &
+                                 values = data_surf_1d) )
+       call check( nf90_put_var( ncid   = ncout_surf_grpid,       &
+                                 varid  = ncout_surf_varids(ivar),   &
+                                 start  = [1],                    & 
+                                 count  = [nsnap],                &
+                                 values = data_surf_1d))
+    end do
     deallocate(data_surf_1d)
 
 

SOLVER/UTILS/post_processing.F90

@@ -282,7 +282,7 @@ program post_processing_seis
   ! define time series
   allocate(time(nt_seis), seis(nt_seis,3))
   do iseis=1, nt_seis
-     time(iseis) = iseis * dt_seis
+     time(iseis) = (iseis - 1) * dt_seis
   enddo
   allocate(seis_sglcomp(nt_seis,3))
 

SOLVER/data_io.f90

@@ -34,19 +34,16 @@ module data_io
   logical           :: dump_snaps_glob
   logical           :: dump_vtk 
   logical           :: dump_xdmf
-  logical           :: dump_snaps_solflu
-  !> N.B. This is not wavefield snapshots, but kernel wavefields. Belongs to nstrain and
-  !! istrain
   logical           :: dump_wavefields 
   logical           :: checkpointing
   logical           :: diagfiles !< Write diagnostic files (seismograms at antipodes, 
                                  !! list of surface elements, blabla), default: false
 
   logical           :: need_fluid_displ
   real(kind=dp)     :: strain_samp
-  integer           :: iseismo  !< current seismogram sample
-  integer           :: istrain  !< current kernel wavefield sample
-  integer           :: isnap    !< current wavefield sample for movies
+  !integer           :: iseismo  !< current seismogram sample
+  !integer           :: istrain  !< current kernel wavefield sample
+  !integer           :: isnap    !< current wavefield sample for movies
   integer           :: nseismo  !< Number of seismogram samples
   integer           :: nstrain  !< Number of wavefield dumps for kernels
   integer           :: nsnap    !< Number of wavefield snapshots for movies

SOLVER/meshes_io.F90

@@ -40,8 +40,8 @@ module meshes_io
   public :: dump_wavefields_mesh_1d
   public :: dump_glob_grid_midpoint
   public :: dump_xdmf_grid
-  public :: dump_solid_grid
-  public :: dump_fluid_grid
+  !public :: dump_solid_grid
+  !public :: dump_fluid_grid
   public :: prepare_mesh_memoryvar_vtk
   public :: build_kwf_grid
   public :: dump_kwf_midpoint_xdmf
@@ -788,14 +788,21 @@ subroutine dump_kwf_grid()
 
   points_mp = 0.
 
+  ct = 1
   do iel=1, nel_solid
-     points_mp(iel,1) = scoord(npol/2,npol/2,ielsolid(iel))
-     points_mp(iel,2) = zcoord(npol/2,npol/2,ielsolid(iel))
+     if (kwf_mask(npol/2,jpol/2,iel)) then
+        points_mp(ct,1) = scoord(npol/2,npol/2,ielsolid(iel))
+        points_mp(ct,2) = zcoord(npol/2,npol/2,ielsolid(iel))
+        ct = ct + 1
+     endif
   enddo
 
   do iel=1, nel_fluid
-     points_mp(iel + nel_solid,1) = scoord(npol/2,npol/2,ielfluid(iel))
-     points_mp(iel + nel_solid,2) = zcoord(npol/2,npol/2,ielfluid(iel))
+     if (kwf_mask(npol/2,jpol/2,iel + nel_solid)) then
+        points_mp(ct,1) = scoord(npol/2,npol/2,ielfluid(iel))
+        points_mp(ct,2) = zcoord(npol/2,npol/2,ielfluid(iel))
+        ct = ct + 1
+     endif
   enddo
 
   if (use_netcdf) then
@@ -1012,25 +1019,25 @@ subroutine dump_kwf_sem_xdmf(filename, npoints, nelem)
 !> Dumps the mesh (s,z) [m] in ASCII format as needed to visualize snapshots 
 !! in the solid region only.
 !! Convention for order in the file: First the fluid, then the solid domain.
-subroutine dump_solid_grid(ibeg,iend,jbeg,jend)
-
-  
-  integer, intent(in) :: ibeg,iend,jbeg,jend 
-  integer             :: iel, ipol,jpol
-
-  open(unit=2500+mynum,file=datapath(1:lfdata)//'/solid_grid_'&
-                            //appmynum//'.dat')
-  do iel=1,nel_solid
-     do jpol=jbeg,jend
-        do ipol=ibeg,iend
-           write(2500+mynum,*)scoord(ipol,jpol,ielsolid(iel)), &
-                              zcoord(ipol,jpol,ielsolid(iel))
-        enddo
-     enddo
-  enddo
-  close(2500+mynum)
-
-end subroutine dump_solid_grid
+!subroutine dump_solid_grid(ibeg,iend,jbeg,jend)
+!
+!  
+!  integer, intent(in) :: ibeg,iend,jbeg,jend 
+!  integer             :: iel, ipol,jpol
+!
+!  open(unit=2500+mynum,file=datapath(1:lfdata)//'/solid_grid_'&
+!                            //appmynum//'.dat')
+!  do iel=1,nel_solid
+!     do jpol=jbeg,jend
+!        do ipol=ibeg,iend
+!           write(2500+mynum,*)scoord(ipol,jpol,ielsolid(iel)), &
+!                              zcoord(ipol,jpol,ielsolid(iel))
+!        enddo
+!     enddo
+!  enddo
+!  close(2500+mynum)
+!
+!end subroutine dump_solid_grid
 !-----------------------------------------------------------------------------------------
 
 !-----------------------------------------------------------------------------------------
@@ -1040,43 +1047,43 @@ end subroutine dump_solid_grid
 !! When reading the fluid wavefield, one therefore needs to multiply all 
 !! components with inv_rho_fluid and the phi component with one/scoord!
 !! Convention for order in the file: First the fluid, then the solid domain.
-subroutine dump_fluid_grid(ibeg,iend,jbeg,jend)
-
-  use data_pointwise, only : inv_rho_fluid
-  
-  
-  integer, intent(in) :: ibeg,iend,jbeg,jend
-  integer             :: iel, ipol,jpol
-  
-  ! When reading the fluid wavefield, one needs to multiply all components 
-  ! with inv_rho_fluid and the phi component with one/scoord!!
-
-  open(unit=2500+mynum,file=datapath(1:lfdata)//&
-                            '/fluid_grid_'//appmynum//'.dat')
-  open(unit=2600+mynum,file=datapath(1:lfdata)//&
-                            '/inv_rho_scoord_fluid_flusnaps_'&
-                            //appmynum//'.dat', STATUS="REPLACE")
-  do iel=1,nel_fluid
-     do jpol=jbeg,jend
-        do ipol=ibeg,iend
-           write(2500+mynum,*)scoord(ipol,jpol,ielfluid(iel)), &
-                              zcoord(ipol,jpol,ielfluid(iel))
-           if ( axis_fluid(iel) .and. ipol==0 ) then
-              ! Axis s=0! write 1 instead of 1/s and then multiply 
-              ! with the correct factor dsdchi, obtained by L'Hospital's rule 
-              ! (see routine fluid_snapshot below).
-              write(2600+mynum,*)inv_rho_fluid(ipol,jpol,iel),one
-           else  
-              write(2600+mynum,*)inv_rho_fluid(ipol,jpol,iel), &
-                                 one/scoord(ipol,jpol,ielfluid(iel))
-           endif
-        enddo
-     enddo
-  enddo
-  close(2500+mynum)
-  close(2600+mynum)
-
-end subroutine dump_fluid_grid
+!subroutine dump_fluid_grid(ibeg,iend,jbeg,jend)
+!
+!  use data_pointwise, only : inv_rho_fluid
+!  
+!  
+!  integer, intent(in) :: ibeg,iend,jbeg,jend
+!  integer             :: iel, ipol,jpol
+!  
+!  ! When reading the fluid wavefield, one needs to multiply all components 
+!  ! with inv_rho_fluid and the phi component with one/scoord!!
+!
+!  open(unit=2500+mynum,file=datapath(1:lfdata)//&
+!                            '/fluid_grid_'//appmynum//'.dat')
+!  open(unit=2600+mynum,file=datapath(1:lfdata)//&
+!                            '/inv_rho_scoord_fluid_flusnaps_'&
+!                            //appmynum//'.dat', STATUS="REPLACE")
+!  do iel=1,nel_fluid
+!     do jpol=jbeg,jend
+!        do ipol=ibeg,iend
+!           write(2500+mynum,*)scoord(ipol,jpol,ielfluid(iel)), &
+!                              zcoord(ipol,jpol,ielfluid(iel))
+!           if ( axis_fluid(iel) .and. ipol==0 ) then
+!              ! Axis s=0! write 1 instead of 1/s and then multiply 
+!              ! with the correct factor dsdchi, obtained by L'Hospital's rule 
+!              ! (see routine fluid_snapshot below).
+!              write(2600+mynum,*)inv_rho_fluid(ipol,jpol,iel),one
+!           else  
+!              write(2600+mynum,*)inv_rho_fluid(ipol,jpol,iel), &
+!                                 one/scoord(ipol,jpol,ielfluid(iel))
+!           endif
+!        enddo
+!     enddo
+!  enddo
+!  close(2500+mynum)
+!  close(2600+mynum)
+!
+!end subroutine dump_fluid_grid
 !-----------------------------------------------------------------------------------------
 
 !-----------------------------------------------------------------------------------------