[cig-commits] r8551 - seismo/2D/SPECFEM2D/trunk

[walter at geodynamics.org]

Author: walter
Date: 2007-12-07 15:55:44 -0800 (Fri, 07 Dec 2007)
New Revision: 8551

Modified:
   seismo/2D/SPECFEM2D/trunk/create_color_image.f90
   seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
added P velocity model as background on color images


Modified: seismo/2D/SPECFEM2D/trunk/create_color_image.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/create_color_image.f90	2007-07-04 00:59:03 UTC (rev 8550)
+++ seismo/2D/SPECFEM2D/trunk/create_color_image.f90	2007-12-07 23:55:44 UTC (rev 8551)
@@ -11,7 +11,7 @@
 !
 !========================================================================
 
-  subroutine create_color_image(color_image_2D_data,iglob_image_color_2D,NX,NY,it,cutsnaps)
+  subroutine create_color_image(color_image_2D_data,iglob_image_color_2D,NX,NY,it,cutsnaps,vp_display,npoin)
 
 ! display a given field as a red and blue color image
 
@@ -23,19 +23,20 @@
 
   include "constants.h"
 
-  integer NX,NY,it
+  integer :: NX,NY,it,npoin
 
-  double precision cutsnaps
+  double precision :: cutsnaps
 
   integer, dimension(NX,NY) :: iglob_image_color_2D
 
   double precision, dimension(NX,NY) :: color_image_2D_data
+  double precision, dimension(npoin) :: vp_display
 
-  integer ix,iy,R,G,B,tenthousands,thousands,hundreds,tens,units,remainder,current_rec
+  integer :: ix,iy,R,G,B,tenthousands,thousands,hundreds,tens,units,remainder,current_rec
 
-  double precision amplitude_max,normalized_value
+  double precision :: amplitude_max,normalized_value,vpmin,vpmax,x1
 
-  character(len=100) file_name,system_command
+  character(len=100) :: file_name,system_command
 
 ! create temporary image files in binary PNM P6 format (smaller) or ASCII PNM P3 format (easier to edit)
   logical, parameter :: BINARY_FILE = .true.
@@ -121,6 +122,8 @@
 
 ! compute maximum amplitude
   amplitude_max = maxval(abs(color_image_2D_data))
+  vpmin = minval(vp_display)
+  vpmax = maxval(vp_display)
 
 ! in the PNM format, the image starts in the upper-left corner
   do iy=NY,1,-1
@@ -137,11 +140,29 @@
 ! suppress small amplitudes considered as noise
       else if (abs(color_image_2D_data(ix,iy)) < amplitude_max * cutsnaps) then
 
-! use black background where amplitude is negligible
-          R = 0
-          G = 0
-          B = 0
+! use P velocity model as background where amplitude is negligible
+        if((vpmax-vpmin)/vpmin > 0.02d0) then
+          x1 = (vp_display(iglob_image_color_2D(ix,iy))-vpmin)/(vpmax-vpmin)
+        else
+          x1 = 0.5d0
+        endif
 
+! rescale to avoid very dark gray levels
+        x1 = x1*0.7 + 0.2
+        if(x1 > 1.d0) x1=1.d0
+
+! invert scale: white = vpmin, dark gray = vpmax
+        x1 = 1.d0 - x1
+
+! map to [0,255]
+        x1 = x1 * 255.d0
+
+        R = nint(x1)
+        if(R < 0) R = 0
+        if(R > 255) R = 255
+        G = R
+        B = R
+
       else
 
 ! define normalized image data in [-1:1] and convert to nearest integer

Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90	2007-07-04 00:59:03 UTC (rev 8550)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90	2007-12-07 23:55:44 UTC (rev 8551)
@@ -144,6 +144,8 @@
   double precision, dimension(:), allocatable :: potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
   double precision, dimension(:), allocatable :: rmass_inverse_elastic,rmass_inverse_acoustic,density,displread,velocread,accelread
 
+  double precision, dimension(:), allocatable :: vp_display
+
   double precision, dimension(:,:,:), allocatable :: vpext,vsext,rhoext
   double precision :: previous_vsext
 
@@ -715,8 +717,6 @@
 
   if(nrec < 1) call exit_MPI('need at least one receiver')
 
-
-
 ! receiver information
   allocate(ispec_selected_rec(nrec))
   allocate(st_xval(nrec))
@@ -1586,6 +1586,26 @@
   time_start = 86400.d0*time_values(3) + 3600.d0*time_values(5) + &
                60.d0*time_values(6) + time_values(7) + time_values(8) / 1000.d0
 
+! to display the P-velocity model in background on color images
+  allocate(vp_display(npoin))
+  do ispec = 1,nspec
+! get relaxed elastic parameters of current spectral element
+    rhol = density(kmato(ispec))
+    lambdal_relaxed = elastcoef(1,kmato(ispec))
+    mul_relaxed = elastcoef(2,kmato(ispec))
+    do j = 1,NGLLZ
+      do i = 1,NGLLX
+!--- if external medium, get elastic parameters of current grid point
+          if(assign_external_model) then
+            vp_display(ibool(i,j,ispec)) = vpext(i,j,ispec)
+          else
+            cpsquare = (lambdal_relaxed + 2.d0*mul_relaxed) / rhol
+            vp_display(ibool(i,j,ispec)) = sqrt(cpsquare)
+          endif
+      enddo
+    enddo
+  enddo
+
 ! *********************************************************
 ! ************* MAIN LOOP OVER THE TIME STEPS *************
 ! *********************************************************
@@ -2135,25 +2155,20 @@
 
      end if
   end if
-!  call MPI_BARRIER(MPI_COMM_WORLD, ier)
 
 #endif
 
 
   if ( myrank == 0 ) then
-     call create_color_image(image_color_data,iglob_image_color,NX_IMAGE_color,NZ_IMAGE_color,it,cutsnaps)
-
+     call create_color_image(image_color_data,iglob_image_color,NX_IMAGE_color,NZ_IMAGE_color,it,cutsnaps,vp_display,npoin)
      write(IOUT,*) 'Color image created'
+  endif
 
-  end if
-
   endif
 
 !----  save temporary or final seismograms
-  if ( it == NSTEP ) then
-     call write_seismograms(sisux,sisuz,station_name,network_name,NSTEP, &
-          nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,it,t0)
-  end if
+  call write_seismograms(sisux,sisuz,station_name,network_name,NSTEP, &
+        nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,it,t0)
 
 ! count elapsed wall-clock time
   call date_and_time(datein,timein,zone,time_values)