SLiMv2.1: SLiMv2.1/aaron_wave.f90 Source File

00001 subroutine wave()
00002   !copyright Robert Cameron. 
00003   include 'comments.f90'
00004   implicit none
00005   include 'param.f90'
00006   include 'terms3d.f90'
00007   include 'functionlist.f90'
00008 
00009   real  t,PE,KE,corr,corr2,FSmax,FSmax1,dt1,random_jitter(nx,ny),rn
00010   complex fdold(nx/2+1,ny,nz_max,nvar)
00011   integer  i,ol,ol_old,i1,k,n,mark,j,rs,rseed(4),nsamp,ssamp    
00012   complex fdtemp(nx/2+1,ny,nz_max,nvar)
00013   complex temp(nx/2+1,ny,nz_max)
00014   integer ii,jj
00015   complex src
00016   real src_r,src_i
00017   integer N_kick
00018   complex the_kick, evolve_a,evolve_b
00019   integer DO_SOURCE
00020   real kick_amp
00021   character*80 src_real_file
00022   character*80 src_imag_file
00023 
00024   parameter (DO_SOURCE=1)
00025 
00026   fsmax=0       
00027   CALL SETUPcoeff(coeff)
00028         do k=1,nz
00029             if(kk(k) .lt. nz_ghostglobal-35) then
00030                FSmax1=sqrt(maxval(&
00031                     & (    coeff(:,:,k,15,1)**2+coeff(:,:,k,16,1)**2&
00032                     &                                +coeff(:,:,k,17,1)**2)&
00033                     &                /(4. *3.1415926)*coeff(:,:,k,12,1)&
00034                     & +coeff(:,:,k,5,1) &
00035                     & ))
00036                !write(*,*),k,kk(k),minval(coeff(:,:,k,5,1)),fsmax1,fsmax
00037             endif
00038             FSmax=max(FSmax,FSmax1)
00039          enddo
00040          dt1=min(dx,dy, dz)/FSmax/4.
00041          call mpi_reduce (dt1, dt, 1, MPI_REAL, MPI_MIN, 0, MPI_COMM_WORLD, ierr )
00042          dt=15./(int(15./dt)+1) 
00043 
00044          call mpi_bcast (dt, 1, MPI_REAL, 0, MPI_COMM_WORLD, ierr)
00045          !full cube sample
00046          nsamp=(30+dt/30.)/(2*dt) 
00047          !slice sample
00048          ssamp=10000*(60+dt/30.)/(2*dt)
00049 
00050          !
00051          !
00052          !source kick time
00053          ! ten seconds is good
00054          ! acb
00055          N_kick=(10+dt/10.)/(2*dt)
00056 
00057          !     INITIALIZE FOURIER TRANSFORMS and derivative arrays and data
00058          call initialize(ider,diff,dt)
00059          call initfftw()
00060          !write(*,*) my_rank,"dt=",dt
00061          CALL input(ol_old)
00062          !ol_old=ol_old*nsamp
00063          ol_old=0
00064          t=ol_old*2*dt
00065          !if (my_rank .eq. 0) write(*,*) "restoring from step:",ol_old, "t=", t
00066   !     -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
00067   step=1
00068 
00069   !
00070   ! start the source work
00071   !
00072   ! acb
00073   write(*,*) "reading source amplitude file"
00074   src_amp_file_name="source_amplitude.fits"
00075   call fr(src_amp_file_name,src_amp_space,nx,ny,nz_ghostglobal+2)
00076   !
00077   ! fourier transform
00078   ! acb
00079   write(*,*) "taking Fourier transform of source amplitude"
00080   call sfftw_execute(plan_src_amp_to_fourier)
00081   !  src_amp_fourier=abs(src_amp_fourier)
00082 
00083 
00084 
00085 
00086   call derivs(PE,KE)
00087   mark=-1
00088   ! the time integration LAX-Wendroff 
00089   do ol=ol_old,10000000   
00090      fdold=fd
00091      step=1
00092      call MPI_ghost()
00093      call BCs()
00094      dfddt(:,:,:,1)=centref(fd(:,:,:,4))
00095      dfddt(:,:,:,2)=centref(fd(:,:,:,5))
00096      dfddt(:,:,:,3)=centref(fd(:,:,:,6))
00097      call derivs(PE,KE)
00098 
00099      if (mod(ol,nsamp) .eq. 0) then
00100         mark=ol/nsamp
00101         if(my_rank .eq. 0) write(*,*) "saving cube",t,mark
00102         call outputcube(fd,rd,rpress,rrhoe,rdb,t,mark)
00103      endif
00104      if (mod(ol,ssamp) .eq. 0) then
00105         mark=ol/ssamp
00106         if (my_rank .eq. 0) write(*,*) "saving slice",t,mark,ke,pe
00107         call outputslice(fd,rd,rpress,rrhoe,rdb,t,mark)
00108      endif
00109      !if(my_rank .eq. 0) write(*,*) t,PE,KE,ol,nsamp
00110 !END AVERAGE
00111      fd(:,:,1:nz-1,:)=(fd(:,:,1:nz-1,:)+fd(:,:,2:nz  ,:))/2.  
00112      fd=fd+(dfddt)*dt  
00113 
00114      ! MPI PASS STUFF HERE
00115      call MPI_ghost()
00116      call BCs() 
00117      step=2
00118      call derivs(PE,KE)
00119      dfddt(:,:,:,1)=(fdold(:,:,:,4))
00120      dfddt(:,:,:,2)=(fdold(:,:,:,5))
00121      dfddt(:,:,:,3)=(fdold(:,:,:,6))
00122 
00123 
00124      !
00125      ! ok, here is the main part of the calculation
00126      ! notice the factor two in front of the dt everywhere  
00127      ! 
00128      ! acb
00129      if (DO_SOURCE .eq. 1) then
00130 
00131         kick_amp=sqrt(2*N_kick*dt / tau)
00132         evolve_a=cmplx(cos(2*coef_b*dt),   sin(2*coef_b*dt))*exp(-2*dt/tau)
00133         evolve_b=cmplx(cos(2*coef_b*dt),-1*sin(2*coef_b*dt))*exp(-2*dt/tau)
00134 
00135         if (mod(ol,N_kick) .eq. 0) then
00136            !
00137            ! if source kick time
00138            !
00139            do ii=1,nx/2+1
00140               do jj=1,ny
00141                  !
00142                  ! get random gaussian numbers
00143                  !
00144                  call GRNF(src_r, src_i)
00145                  the_kick=kick_amp*cmplx(src_r,src_i)
00146                  src_fourier_a(ii,jj,:)=src_fourier_a(ii,jj,:)+the_kick
00147                  
00148                  call GRNF(src_r, src_i)           
00149                  the_kick=kick_amp*cmplx(src_r,src_i)
00150                  src_fourier_b(ii,jj,:)=src_fourier_b(ii,jj,:)+the_kick
00151              
00152               end do
00153            end do
00154         else
00155            !
00156            ! else just evolve the source fcns
00157            !
00158            ! write(*,*) "Evolve", t, evolve
00159            do ii=1,nx/2+1
00160               do jj=1,ny
00161                  src_fourier_a(ii,jj,:)=src_fourier_a(ii,jj,:)*evolve_a
00162                  src_fourier_b(ii,jj,:)=src_fourier_b(ii,jj,:)*evolve_b
00163               end do
00164            end do
00165         end if
00166         
00167     
00168         
00169 
00170         !
00171         ! now add the forcing function to the vertical displacement
00172         ! acb
00173         !
00174         dfddt(:,:,:,3)= dfddt(:,:,:,3)+src_fourier_a*src_amp_fourier
00175         dfddt(:,:,:,3)= dfddt(:,:,:,3)+src_fourier_b*src_amp_fourier
00176 
00177         
00178      endif
00179      
00180      
00181 
00182 
00183 
00184 
00185      fd=fdold+(dfddt)*2.*dt
00186      t=t+2.*dt
00187      call MPI_ghost()
00188      !===========================
00189      do i1=1,6
00190           temp=d2fdz2(fd(:,:,:,i1))
00191           diffdfddt(:,:,:,i1)=(temp)
00192      enddo
00193      fd(:,:,:,:)=fd(:,:,:,:)+diffdfddt(:,:,:,:)*2*dt*nu*2000   
00194      do i1=1,6
00195         !
00196         ! I stuck this multiple by 4 here to increase the
00197         ! damping the boundary layers
00198         !
00199         ! feb 26, 2008
00200         !
00201         ! acb
00202              diffdfddt(:,:,:,i1)=fd(:,:,:,i1)*diff(:,:,:)*4
00203      enddo
00204      fd=fd+diffdfddt*2*dt
00205      call mpi_barrier(MPI_COMM_WORLD,ierr)
00206   enddo
00207   !     =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
00208 end subroutine wave
00209 
00210 !================================================================
00211 
00212 subroutine BCs()
00213   include 'comments.f90'
00214   implicit none
00215   include 'param.f90'
00216   include 'terms3d.f90'
00217   include 'functionlist.f90'
00218   complex  tempin(nx/2+1,ny,2,nvar)
00219   complex tempout(nx/2+1,ny,2,nvar)
00220   integer i     
00221   if(isbottom) then
00222      do i=1,6 
00223         fd(:,:,1,i)=fd(:,:,2,i)*0.99985
00224     enddo
00225   endif   
00226   if(istop) then
00227      fd(:,:,nz,1:6)=(2.*fd(:,:,nz-1,1:6)-fd(:,:,nz-2,1:6))
00228      fd(:,:,nz-1,1:6)=(2.*fd(:,:,nz-1,1:6)-fd(:,:,nz-2,1:6))
00229 !     fd(:,:,nz-20:nz,1:6)=0
00230   endif
00231 end subroutine BCs
00232 !============================================
00233 subroutine MPI_ghost()
00234   include 'comments.f90'
00235   implicit none
00236   include 'param.f90'
00237   include 'terms3d.f90'
00238   include 'functionlist.f90'
00239   complex  tempin(nx/2+1,ny,2,nvar)
00240   complex tempout(nx/2+1,ny,2,nvar)
00241   integer  nmess,up,down,ms
00242   integer*8 i
00243   ms=(nx/2+1)*ny*2*nvar
00244   up=my_rank+1
00245   down=my_rank-1
00246   nmess=12
00247   tempout=fd(:,:,nz-3:nz-2,:)
00248   if(istop .eqv. .false.) then
00249      call mpi_send(tempout,ms,MPI_COMPLEX,&
00250           &                         up,nmess,MPI_COMM_WORLD,ierr)
00251   endif
00252   if(isbottom .eqv. .false.) then
00253              call mpi_recv(tempin ,ms,MPI_COMPLEX,&
00254           &                         down,nmess,MPI_COMM_WORLD,istatus,ierr)
00255              fd(:,:,1:2,:)=tempin
00256   endif 
00257 
00258 
00259   tempout=fd(:,:,3:4,:)
00260   if(isbottom .eqv. .false.) call mpi_send(tempout,ms,MPI_COMPLEX,&
00261           &                         down,nmess,MPI_COMM_WORLD,ierr)
00262   if(istop .eqv. .false.)   then 
00263      call mpi_recv(tempin ,ms,MPI_COMPLEX,&
00264           &                         up,nmess,MPI_COMM_WORLD,istatus,ierr)
00265      fd(:,:,nz-1:nz,:)=tempin
00266   endif  
00267 !  call mpi_finalize(ierr)
00268 !  stop
00269 
00270 
00271 end subroutine MPI_ghost
00272 
00273 
00274 subroutine fr(fn,data,nx,ny,nz)
00275   implicit none 
00276   integer nx,ny,nz,blocksize    
00277   character*80 fn
00278   real data(nx,ny,nz)
00279   integer readwrite,group,fpixel,nelements,status,unit
00280   logical er
00281   character*30 errmsg   
00282 
00283   readwrite=0   
00284   blocksize=1
00285   status=0      
00286   call ftgiou(unit,status)      
00287   call ftopen(unit,fn,readwrite,blocksize,status)
00288   group=1
00289   fpixel=1
00290   nelements=nx*ny*nz
00291   call ftgpve(unit,group,fpixel,nelements,0.,data,er,status)
00292   call ftclos(unit, status)
00293   call ftfiou(unit, status)
00294 
00295   if(status .ne. 0) then
00296      call ftgerr(status,errmsg)
00297      write(*,*) "reading: ",fn
00298      write(*,*) "readin error: ",status,errmsg
00299      stop
00300   endif
00301 end subroutine fr
00302 
00303 subroutine fw(fn,data,nx,ny,nz,n_rec)
00304   implicit none 
00305   character*80 fn,record,infn
00306   character*30 errmsg
00307   real data(nx,ny,nz),alpha,b_0 
00308   integer status,unit,blocksize,bitpix,naxis,naxes(3)
00309   integer i,j,group,fpixel,nelements,nx,ny,nz
00310   integer n,n_rec,readwrite,nkeys,nspace,iunit
00311   logical simple,extend
00312   
00313   status=0
00314   write(*,*) n_rec,minval(data),maxval(data)    
00315   !     Get an unused Logical Unit Number to use to create the FITS file
00316   call ftgiou(unit,status)
00317   call  ftgiou(iunit,status)
00318   blocksize=1
00319   readwrite=1  
00320   if(n_rec .eq. 0)      then
00321      call ftinit(unit,fn,blocksize,status)
00322   else 
00323      call ftnopn(unit,fn,blocksize,status)
00324      call ftcrhd(unit,status)
00325   endif
00326   simple=.true.
00327   bitpix=-32
00328   naxis=3
00329   naxes(1)=nx
00330   naxes(2)=ny
00331   naxes(3)=nz   
00332   extend=.true.
00333   call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00334   
00335   if (status .ne. 0) then
00336      call ftgerr(status,errmsg) 
00337      write(*,*) "error:",errmsg
00338      stop
00339   endif
00340   !OK now we need to open a data file to retrieve the header information
00341   call  ftgiou(iunit,status)
00342   readwrite=0
00343   infn= "spot/bx1_wavesim.fits"
00344   call ftopen(iunit,infn,readwrite,blocksize,status)
00345   !how many key words do we have?
00346   call ftghsp(iunit,nkeys,nspace,status)
00347   ! Read each 80-character keyword record, and print it out.
00348   do i = 9, nkeys
00349      call ftgrec(iunit,i,record,status)
00350      call ftprec(unit,record,status)
00351   enddo
00352   call ftclos(iunit, status)
00353   call ftfiou(iunit, status)
00354   !OK HEADER DONE
00355   !     write the array to the FITS file
00356   group=1
00357   fpixel=1
00358   nelements=nx*ny*nz
00359   call ftppre(unit,group,fpixel,nelements,data,status)  
00360   !make next hdu
00361   call ftclos(unit, status)
00362   call ftfiou(unit, status)     
00363 end subroutine fw
00364 
00365 
00366 SUBROUTINE GRNF (X,Y)
00367 !
00368 ! Two Gaussian random numbers generated from two uniform random
00369 ! numbers. Copyright (c) Tao Pang 1997.
00370 !
00371   IMPLICIT NONE
00372   REAL, INTENT (OUT) :: X,Y
00373   REAL :: PI,R1,R2,R,RANF
00374 !
00375   PI = 4.0*ATAN(1.0)
00376   R1 = -ALOG(1.0-RAND())
00377   R2 = 2.0*PI*RAND()
00378   R1 = SQRT(2.0*R1)
00379   X  = R1*COS(R2)
00380   Y  = R1*SIN(R2)
00381 END SUBROUTINE GRNF