[cig-commits] r18171 - in seismo/3D/SPECFEM3D_GLOBE/trunk/src: shared specfem3D

danielpeter at geodynamics.org danielpeter at geodynamics.org
Mon Apr 4 14:33:00 PDT 2011


Author: danielpeter
Date: 2011-04-04 14:33:00 -0700 (Mon, 04 Apr 2011)
New Revision: 18171

Modified:
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_c_binary.c
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90
Log:
updates noise routines to take pre-allocated surface movie array

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_c_binary.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_c_binary.c	2011-04-04 18:28:44 UTC (rev 18170)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_c_binary.c	2011-04-04 21:33:00 UTC (rev 18171)
@@ -289,6 +289,7 @@
     if (ret > 0){
       donelen = donelen + ret;
       remlen = remlen - MAX_B;
+      // this shifts the pointer position, thus is used in arithmetic, compilers might warn about this...
       buf += MAX_B;
     }
     else{
@@ -334,6 +335,7 @@
     if (ret > 0){
       donelen = donelen + ret;
       remlen = remlen - MAX_B;
+      // this shifts the pointer position, thus is used in arithmetic, compilers might warn about this...
       buf += MAX_B;
     }
     else{

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90	2011-04-04 18:28:44 UTC (rev 18170)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90	2011-04-04 21:33:00 UTC (rev 18171)
@@ -82,10 +82,10 @@
                                    xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
                                    irec_master_noise,normal_x_noise,normal_y_noise,normal_z_noise,mask_noise)
   implicit none
+  include 'mpif.h'
+  include "precision.h"
   include "constants.h"
   include "OUTPUT_FILES/values_from_mesher.h"
-  include 'mpif.h'
-  include "precision.h"
   ! input parameters
   integer :: myrank, nrec, NSTEP, nmovie_points, nspec_top, NIT
   integer, dimension(nrec) :: islice_selected_rec
@@ -121,9 +121,12 @@
   close(IIN_NOISE)
 
   if (myrank == 0) then
-     open(unit=IOUT_NOISE,file='OUTPUT_FILES/irec_master_noise',status='unknown',action='write')
-     WRITE(IOUT_NOISE,*) 'The master receiver is: (RECEIVER ID)', irec_master_noise
-     close(IOUT_NOISE)
+    open(unit=IOUT_NOISE,file='OUTPUT_FILES/irec_master_noise', &
+          status='unknown',action='write',iostat=ios)
+    if( ios /= 0 ) call exit_MPI(myrank,'error opening output file irec_master_noise')
+     
+    WRITE(IOUT_NOISE,*) 'The master receiver is: (RECEIVER ID)', irec_master_noise
+    close(IOUT_NOISE)
   endif
 
   ! compute source arrays for "ensemble forward source", which is source of "ensemble forward wavefield"
@@ -160,48 +163,52 @@
   enddo
 
   !!!BEGIN!!! save mask_noise for check, a file called "mask_noise" is saved in "./OUTPUT_FIELS/"
-    ipoin = 0
-      do ispec2D = 1, nspec_top ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
-          ispec = ibelm_top_crust_mantle(ispec2D)
-          k = NGLLZ
-        ! loop on all the points inside the element
-          do j = 1,NGLLY,NIT
-             do i = 1,NGLLX,NIT
-                ipoin = ipoin + 1
-                iglob = ibool_crust_mantle(i,j,k,ispec)
-                store_val_x(ipoin) = xstore_crust_mantle(iglob)
-                store_val_y(ipoin) = ystore_crust_mantle(iglob)
-                store_val_z(ipoin) = zstore_crust_mantle(iglob)
-                store_val_ux(ipoin) = mask_noise(ipoin)
-                store_val_uy(ipoin) = mask_noise(ipoin)
-                store_val_uz(ipoin) = mask_noise(ipoin)
-             enddo
-          enddo
+  ipoin = 0
+  do ispec2D = 1, nspec_top ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+    ispec = ibelm_top_crust_mantle(ispec2D)
+    k = NGLLZ
+    ! loop on all the points inside the element
+    do j = 1,NGLLY,NIT
+      do i = 1,NGLLX,NIT
+        ipoin = ipoin + 1
+        iglob = ibool_crust_mantle(i,j,k,ispec)
+        store_val_x(ipoin) = xstore_crust_mantle(iglob)
+        store_val_y(ipoin) = ystore_crust_mantle(iglob)
+        store_val_z(ipoin) = zstore_crust_mantle(iglob)
+        store_val_ux(ipoin) = mask_noise(ipoin)
+        store_val_uy(ipoin) = mask_noise(ipoin)
+        store_val_uz(ipoin) = mask_noise(ipoin)
       enddo
+    enddo
+  enddo
 
   ! gather info on master proc
-      ispec = nmovie_points
-      call MPI_GATHER(store_val_x,ispec,CUSTOM_MPI_TYPE,store_val_x_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
-      call MPI_GATHER(store_val_y,ispec,CUSTOM_MPI_TYPE,store_val_y_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
-      call MPI_GATHER(store_val_z,ispec,CUSTOM_MPI_TYPE,store_val_z_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
-      call MPI_GATHER(store_val_ux,ispec,CUSTOM_MPI_TYPE,store_val_ux_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
-      call MPI_GATHER(store_val_uy,ispec,CUSTOM_MPI_TYPE,store_val_uy_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
-      call MPI_GATHER(store_val_uz,ispec,CUSTOM_MPI_TYPE,store_val_uz_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+  ispec = nmovie_points
+  call MPI_GATHER(store_val_x,ispec,CUSTOM_MPI_TYPE,store_val_x_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+  call MPI_GATHER(store_val_y,ispec,CUSTOM_MPI_TYPE,store_val_y_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+  call MPI_GATHER(store_val_z,ispec,CUSTOM_MPI_TYPE,store_val_z_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+  call MPI_GATHER(store_val_ux,ispec,CUSTOM_MPI_TYPE,store_val_ux_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+  call MPI_GATHER(store_val_uy,ispec,CUSTOM_MPI_TYPE,store_val_uy_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+  call MPI_GATHER(store_val_uz,ispec,CUSTOM_MPI_TYPE,store_val_uz_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
 
   ! save maks_noise data to disk in home directory
   ! this file can be viewed the same way as surface movie data (xcreate_movie_AVS_DX)
   ! create_movie_AVS_DX.f90 needs to be modified in order to do that,
   ! i.e., instead of showing the normal component, change it to either x, y or z component, or the norm.
-    if(myrank == 0) then
-        open(unit=IOUT_NOISE,file='OUTPUT_FILES/mask_noise',status='unknown',form='unformatted',action='write')
-        write(IOUT_NOISE) store_val_x_all
-        write(IOUT_NOISE) store_val_y_all
-        write(IOUT_NOISE) store_val_z_all
-        write(IOUT_NOISE) store_val_ux_all
-        write(IOUT_NOISE) store_val_uy_all
-        write(IOUT_NOISE) store_val_uz_all
-        close(IOUT_NOISE)
-     endif
+  if(myrank == 0) then
+    open(unit=IOUT_NOISE,file='OUTPUT_FILES/mask_noise', &
+              status='unknown',form='unformatted',action='write',iostat=ios)
+    if( ios /= 0 ) call exit_MPI(myrank,'error opening output file mask_noise')
+        
+    write(IOUT_NOISE) store_val_x_all
+    write(IOUT_NOISE) store_val_y_all
+    write(IOUT_NOISE) store_val_z_all
+    write(IOUT_NOISE) store_val_ux_all
+    write(IOUT_NOISE) store_val_uy_all
+    write(IOUT_NOISE) store_val_uz_all
+    
+    close(IOUT_NOISE)
+  endif
   !!!END!!! save mask_noise for check, a file called "mask_noise" is saved in "./OUTPUT_FIELS/"
 
   end subroutine read_parameters_noise
@@ -217,22 +224,26 @@
                                     SAVE_ALL_SEISMOS_IN_ONE_FILE, USE_BINARY_FOR_LARGE_FILE, &
                                     MOVIE_COARSE,LOCAL_PATH,NSPEC_TOP,NSTEP)
   implicit none
-  include 'mpif.h'
-  include "precision.h"
+!  include 'mpif.h'
+!  include "precision.h"
   include "constants.h"
   include "OUTPUT_FILES/values_from_mesher.h"
   ! input parameters
   integer :: myrank,NOISE_TOMOGRAPHY,SIMULATION_TYPE,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NSPEC_TOP,NSTEP
-  logical :: SAVE_FORWARD,ROTATE_SEISMOGRAMS_RT,SAVE_ALL_SEISMOS_IN_ONE_FILE, USE_BINARY_FOR_LARGE_FILE,MOVIE_COARSE
+  logical :: SAVE_FORWARD,ROTATE_SEISMOGRAMS_RT,SAVE_ALL_SEISMOS_IN_ONE_FILE, &
+            USE_BINARY_FOR_LARGE_FILE,MOVIE_COARSE
   character(len=150) :: LOCAL_PATH
   ! output parameters
   ! local parameters
-  integer :: reclen
+  integer :: reclen,ier
   character(len=150) :: outputname
 
 
   if (myrank == 0) then
-     open(unit=IOUT_NOISE,file='OUTPUT_FILES/NOISE_SIMULATION',status='unknown',action='write')
+     open(unit=IOUT_NOISE,file='OUTPUT_FILES/NOISE_SIMULATION', &
+          status='unknown',action='write',iostat=ier)
+     if( ier /= 0 ) call exit_MPI(myrank,'error opening output file NOISE_SIMULATION')
+     
      WRITE(IOUT_NOISE,*) '*******************************************************************************'
      WRITE(IOUT_NOISE,*) '*******************************************************************************'
      WRITE(IOUT_NOISE,*) 'WARNING!!!!!!!!!!!!'
@@ -506,11 +517,12 @@
                     xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
                     store_val_x,store_val_y,store_val_z, &
                     store_val_ux,store_val_uy,store_val_uz, &
-                    ibelm_top_crust_mantle,ibool_crust_mantle,nspec_top, &
+                    ibelm_top_crust_mantle,ibool_crust_mantle, &
+                    nspec_top,noise_surface_movie, &
                     NIT,it,LOCAL_PATH)
   implicit none
-  include 'mpif.h'
-  include "precision.h"
+!  include 'mpif.h'
+!  include "precision.h"
   include "constants.h"
   include "OUTPUT_FILES/values_from_mesher.h"
   ! input parameters
@@ -527,14 +539,13 @@
   real(kind=CUSTOM_REAL), dimension(nmovie_points) :: &
       store_val_x,store_val_y,store_val_z, &
       store_val_ux,store_val_uy,store_val_uz
-  real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: SURFACE_MOVIE
+  real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,nspec_top) :: noise_surface_movie
+  
   !integer :: ipoin
   !character(len=150) :: outputname
   integer :: idummy
   real(kind=CUSTOM_REAL) :: rdummy
   character :: cdummy    
-  
-  allocate(SURFACE_MOVIE(NDIM,NGLLX,NGLLY,nspec_top))
 
   ! get coordinates of surface mesh and surface displacement
   do ispec2D = 1, nspec_top ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
@@ -543,7 +554,7 @@
     do j = 1,NGLLY,NIT
       do i = 1,NGLLX,NIT
         iglob = ibool_crust_mantle(i,j,k,ispec)
-        SURFACE_MOVIE(:,i,j,ispec2D)=displ_crust_mantle(:,iglob)
+        noise_surface_movie(:,i,j,ispec2D) = displ_crust_mantle(:,iglob)
       enddo
     enddo
   enddo
@@ -552,8 +563,7 @@
   ! LOCAL storage is better than GLOBAL, because we have to save the 'movie' at every time step
   ! also note that the surface movie does NOT have to be shared with other nodes/CPUs
   ! change LOCAL_PATH specified in "DATA/Par_file"
-  call write_abs(9,SURFACE_MOVIE,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
-  deallocate(SURFACE_MOVIE)
+  call write_abs(9,noise_surface_movie,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
 
   ! just to avoid compiler warnings
   idummy = myrank
@@ -566,7 +576,7 @@
   rdummy = store_val_ux(1)
   rdummy = store_val_uy(1)
   rdummy = store_val_uz(1)
-  cdummy = LOCAL_PATH(1)
+  cdummy = LOCAL_PATH(1:1)
 
   end subroutine noise_save_surface_movie
 
@@ -659,11 +669,12 @@
   subroutine noise_read_add_surface_movie(myrank,nmovie_points,accel_crust_mantle, &
                   normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, &
                   store_val_ux,store_val_uy,store_val_uz, &
-                  ibelm_top_crust_mantle,ibool_crust_mantle,nspec_top, &
+                  ibelm_top_crust_mantle,ibool_crust_mantle, &
+                  nspec_top,noise_surface_movie, &
                   NIT,it,LOCAL_PATH,jacobian2D_top_crust_mantle,wgllwgll_xy)
   implicit none
-  include 'mpif.h'
-  include "precision.h"
+!  include 'mpif.h'
+!  include "precision.h"
   include "constants.h"
   include "OUTPUT_FILES/values_from_mesher.h"
   ! input parameters
@@ -680,15 +691,16 @@
   integer :: ipoin,ispec2D,ispec,i,j,k,iglob !,ios
   real(kind=CUSTOM_REAL), dimension(nmovie_points) :: store_val_ux,store_val_uy,store_val_uz
   real(kind=CUSTOM_REAL) :: eta
-  real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: SURFACE_MOVIE
+  real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,nspec_top) :: noise_surface_movie
+
   !character(len=150) :: outputname
   integer :: idummy
   real(kind=CUSTOM_REAL) :: rdummy
   character :: cdummy
 
-  allocate(SURFACE_MOVIE(NDIM,NGLLX,NGLLY,nspec_top))
   ! read surface movie
-  call read_abs(9,SURFACE_MOVIE,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
+  call read_abs(9,noise_surface_movie,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
+
   ! get coordinates of surface mesh and surface displacement
   ipoin = 0
   do ispec2D = 1, nspec_top ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
@@ -702,10 +714,11 @@
         ipoin = ipoin + 1
         iglob = ibool_crust_mantle(i,j,k,ispec)
 
-        eta = SURFACE_MOVIE(1,i,j,ispec2D) * normal_x_noise(ipoin) + &
-              SURFACE_MOVIE(2,i,j,ispec2D) * normal_y_noise(ipoin) + &
-              SURFACE_MOVIE(3,i,j,ispec2D) * normal_z_noise(ipoin)
+        eta = noise_surface_movie(1,i,j,ispec2D) * normal_x_noise(ipoin) + &
+              noise_surface_movie(2,i,j,ispec2D) * normal_y_noise(ipoin) + &
+              noise_surface_movie(3,i,j,ispec2D) * normal_z_noise(ipoin)
 
+
         accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + eta * mask_noise(ipoin) * normal_x_noise(ipoin) &
                                                       * wgllwgll_xy(i,j) * jacobian2D_top_crust_mantle(i,j,ispec2D)
         accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + eta * mask_noise(ipoin) * normal_y_noise(ipoin) &
@@ -716,15 +729,13 @@
     enddo
 
   enddo
-  deallocate(SURFACE_MOVIE)
 
-
   ! just to avoid compiler warnings
   idummy = myrank
   rdummy = store_val_ux(1)
   rdummy = store_val_uy(1)
   rdummy = store_val_uz(1)
-  cdummy = LOCAL_PATH(1)
+  cdummy = LOCAL_PATH(1:1)
 
   end subroutine noise_read_add_surface_movie
 
@@ -813,7 +824,8 @@
   subroutine compute_kernels_strength_noise(myrank,ibool_crust_mantle, &
                           Sigma_kl_crust_mantle,displ_crust_mantle,deltat,it, &
                           nmovie_points,normal_x_noise,normal_y_noise,normal_z_noise, &
-                          nspec_top,ibelm_top_crust_mantle,LOCAL_PATH)
+                          nspec_top,noise_surface_movie, &
+                          ibelm_top_crust_mantle,LOCAL_PATH)
   implicit none
   include "constants.h"
   include "OUTPUT_FILES/values_from_mesher.h"
@@ -832,15 +844,16 @@
   integer :: i,j,k,ispec,iglob,ipoin,ispec2D !,ios
   real(kind=CUSTOM_REAL) :: eta
   real(kind=CUSTOM_REAL), dimension(nmovie_points) :: store_val_ux,store_val_uy,store_val_uz
-  real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: SURFACE_MOVIE
+  real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,nspec_top) :: noise_surface_movie
+  
   !character(len=150) :: outputname
   integer :: idummy
   real(kind=CUSTOM_REAL) :: rdummy
   character :: cdummy
   
-  allocate(SURFACE_MOVIE(NDIM,NGLLX,NGLLY,nspec_top))
   ! read surface movie, needed for Sigma_kl_crust_mantle
-  call read_abs(9,SURFACE_MOVIE,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
+  call read_abs(9,noise_surface_movie,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
+
   ! noise source strength kernel
   ! to keep similar structure to other kernels, the source strength kernel is saved as a volumetric kernel
   ! but only updated at the surface, because the noise is generated there
@@ -856,9 +869,9 @@
         ipoin = ipoin + 1
         iglob = ibool_crust_mantle(i,j,k,ispec)
 
-        eta = SURFACE_MOVIE(1,i,j,ispec2D) * normal_x_noise(ipoin) + &
-              SURFACE_MOVIE(2,i,j,ispec2D) * normal_y_noise(ipoin) + &
-              SURFACE_MOVIE(3,i,j,ispec2D) * normal_z_noise(ipoin)
+        eta = noise_surface_movie(1,i,j,ispec2D) * normal_x_noise(ipoin) + &
+              noise_surface_movie(2,i,j,ispec2D) * normal_y_noise(ipoin) + &
+              noise_surface_movie(3,i,j,ispec2D) * normal_z_noise(ipoin)
 
         Sigma_kl_crust_mantle(i,j,k,ispec) =  Sigma_kl_crust_mantle(i,j,k,ispec) &
            + deltat * eta * ( normal_x_noise(ipoin) * displ_crust_mantle(1,iglob) &
@@ -874,9 +887,8 @@
   rdummy = store_val_ux(1)
   rdummy = store_val_uy(1)
   rdummy = store_val_uz(1)
-  cdummy = LOCAL_PATH(1)
+  cdummy = LOCAL_PATH(1:1)
 
-
   end subroutine compute_kernels_strength_noise
 
 ! =============================================================================================================
@@ -896,13 +908,17 @@
   ! output parameters
   ! local parameters
   character(len=150) :: prname
+  integer :: ier
 
-
   call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_PATH)
 
-    open(unit=IOUT_NOISE,file=trim(prname)//'Sigma_kernel.bin',status='unknown',form='unformatted',action='write')
-    write(IOUT_NOISE) Sigma_kl_crust_mantle     ! need to put dimensions back (not done yet)
-    close(IOUT_NOISE)
+  open(unit=IOUT_NOISE,file=trim(prname)//'Sigma_kernel.bin', &
+        status='unknown',form='unformatted',action='write',iostat=ier)
+  if( ier /= 0 ) call exit_MPI(myrank,'error opening file Sigma_kernel.bin')
+  
+  write(IOUT_NOISE) Sigma_kl_crust_mantle     ! need to put dimensions back (not done yet)
+  close(IOUT_NOISE)
+  
   end subroutine save_kernels_strength_noise
 
 ! =============================================================================================================

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90	2011-04-04 18:28:44 UTC (rev 18170)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90	2011-04-04 21:33:00 UTC (rev 18171)
@@ -869,9 +869,10 @@
 
 ! NOISE_TOMOGRAPHY
   real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: noise_sourcearray
-  integer :: irec_master_noise
   real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
              normal_x_noise,normal_y_noise,normal_z_noise, mask_noise
+  real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: noise_surface_movie  
+  integer :: irec_master_noise
 
 ! ************** PROGRAM STARTS HERE **************
 !
@@ -1837,33 +1838,35 @@
                     b_A_array_rotation,b_B_array_rotation,LOCAL_PATH)
 
 !<YANGL
-    ! NOISE TOMOGRAPHY
-    if ( NOISE_TOMOGRAPHY /= 0 ) then
-      allocate(noise_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NSTEP), &
-              normal_x_noise(nmovie_points), &
-              normal_y_noise(nmovie_points), &
-              normal_z_noise(nmovie_points), &
-              mask_noise(nmovie_points),stat=ier)
-       if( ier /= 0 ) call exit_MPI(myrank,'error allocating noise arrays')
+  ! NOISE TOMOGRAPHY
+  if ( NOISE_TOMOGRAPHY /= 0 ) then
+    allocate(noise_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NSTEP), &
+            normal_x_noise(nmovie_points), &
+            normal_y_noise(nmovie_points), &
+            normal_z_noise(nmovie_points), &
+            mask_noise(nmovie_points), &
+            noise_surface_movie(NDIM,NGLLX,NGLLY,NSPEC2D_TOP(IREGION_CRUST_MANTLE)),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating noise arrays')
   
-       noise_sourcearray(:,:,:,:,:) = 0._CUSTOM_REAL
-       normal_x_noise(:)            = 0._CUSTOM_REAL
-       normal_y_noise(:)            = 0._CUSTOM_REAL
-       normal_z_noise(:)            = 0._CUSTOM_REAL
-       mask_noise(:)                = 0._CUSTOM_REAL
+    noise_sourcearray(:,:,:,:,:) = 0._CUSTOM_REAL
+    normal_x_noise(:)            = 0._CUSTOM_REAL
+    normal_y_noise(:)            = 0._CUSTOM_REAL
+    normal_z_noise(:)            = 0._CUSTOM_REAL
+    mask_noise(:)                = 0._CUSTOM_REAL
+    noise_surface_movie(:,:,:,:) = 0._CUSTOM_REAL
+    
+    call read_parameters_noise(myrank,nrec,NSTEP,nmovie_points, &
+                              islice_selected_rec,xi_receiver,eta_receiver,gamma_receiver,nu, &
+                              noise_sourcearray,xigll,yigll,zigll,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
+                              NIT, ibool_crust_mantle, ibelm_top_crust_mantle, &
+                              xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
+                              irec_master_noise,normal_x_noise,normal_y_noise,normal_z_noise,mask_noise)
 
-       call read_parameters_noise(myrank,nrec,NSTEP,nmovie_points, &
-                                  islice_selected_rec,xi_receiver,eta_receiver,gamma_receiver,nu, &
-                                  noise_sourcearray,xigll,yigll,zigll,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
-                                  NIT, ibool_crust_mantle, ibelm_top_crust_mantle, &
-                                  xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
-                                  irec_master_noise,normal_x_noise,normal_y_noise,normal_z_noise,mask_noise)
-
-       call check_parameters_noise(myrank,NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, &
-                                  NUMBER_OF_RUNS, NUMBER_OF_THIS_RUN,ROTATE_SEISMOGRAMS_RT, &
-                                  SAVE_ALL_SEISMOS_IN_ONE_FILE, USE_BINARY_FOR_LARGE_FILE, &
-                                  MOVIE_COARSE,LOCAL_PATH,NSPEC2D_TOP(IREGION_CRUST_MANTLE),NSTEP)
-    endif
+    call check_parameters_noise(myrank,NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, &
+                              NUMBER_OF_RUNS, NUMBER_OF_THIS_RUN,ROTATE_SEISMOGRAMS_RT, &
+                              SAVE_ALL_SEISMOS_IN_ONE_FILE, USE_BINARY_FOR_LARGE_FILE, &
+                              MOVIE_COARSE,LOCAL_PATH,NSPEC2D_TOP(IREGION_CRUST_MANTLE),NSTEP)
+  endif
 !>YANGL
 
 !
@@ -3067,7 +3070,8 @@
        call noise_read_add_surface_movie(myrank,nmovie_points,accel_crust_mantle, &
                               normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, &
                               store_val_ux,store_val_uy,store_val_uz, &
-                              ibelm_top_crust_mantle,ibool_crust_mantle,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
+                              ibelm_top_crust_mantle,ibool_crust_mantle, &
+                              NSPEC2D_TOP(IREGION_CRUST_MANTLE),noise_surface_movie, &
                               NIT,NSTEP-it+1,LOCAL_PATH,jacobian2D_top_crust_mantle,wgllwgll_xy)
         ! be careful, since ensemble forward sources are reversals of generating wavefield "eta"
         ! hence the "NSTEP-it+1", i.e., start reading from the last timestep
@@ -3082,7 +3086,8 @@
         call noise_read_add_surface_movie(myrank,nmovie_points,b_accel_crust_mantle, &
                               normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, &
                               store_val_ux,store_val_uy,store_val_uz, &
-                              ibelm_top_crust_mantle,ibool_crust_mantle,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
+                              ibelm_top_crust_mantle,ibool_crust_mantle, &
+                              NSPEC2D_TOP(IREGION_CRUST_MANTLE),noise_surface_movie, &
                               NIT,it,LOCAL_PATH,jacobian2D_top_crust_mantle,wgllwgll_xy)
     endif
 !>YANGL
@@ -4019,7 +4024,8 @@
        call compute_kernels_strength_noise(myrank,ibool_crust_mantle, &
                           Sigma_kl_crust_mantle,displ_crust_mantle,deltat,it, &
                           nmovie_points,normal_x_noise,normal_y_noise,normal_z_noise, &
-                          NSPEC2D_TOP(IREGION_CRUST_MANTLE),ibelm_top_crust_mantle,LOCAL_PATH)
+                          NSPEC2D_TOP(IREGION_CRUST_MANTLE),noise_surface_movie, &
+                          ibelm_top_crust_mantle,LOCAL_PATH)
 !>YANGL
 
     ! --- boundary kernels ------
@@ -4252,7 +4258,7 @@
                             store_val_x,store_val_y,store_val_z, &
                             store_val_ux,store_val_uy,store_val_uz, &
                             ibelm_top_crust_mantle,ibool_crust_mantle, &
-                            NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
+                            NSPEC2D_TOP(IREGION_CRUST_MANTLE),noise_surface_movie, &
                             NIT,it,LOCAL_PATH)
   endif
 !>YANGL
@@ -4392,7 +4398,10 @@
   endif
 
   ! save/read the surface movie using the same c routine as we do for absorbing boundaries (file ID is 9)
-  if (NOISE_TOMOGRAPHY/=0) call close_file_abs(9) 
+  if (NOISE_TOMOGRAPHY/=0) then
+    call close_file_abs(9) 
+    deallocate(noise_surface_movie)
+  endif
 
 
   ! synchronize all processes



More information about the CIG-COMMITS mailing list