SLiMv2.1: SLiMv2.1/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 
00015   !Aarons stuff
00016   integer ii,jj,ll
00017   real src_r,src_i
00018 
00019   integer N_kick
00020   complex the_kick, evolve_a,evolve_b
00021   integer DO_SOURCE
00022   !need to define these
00023   real tau,t_src
00024   integer n_src
00025   parameter(n_src=60)
00026   real x_src,y_src,tt(n_src),t_corr,told
00027   complex src_move(nx/2+1,ny,n_src),src_fourier(nx/2+1,ny,nz_max)
00028   INTEGER RAN_SEED
00029   !forcing properties
00030   tt=1.e5
00031   src_move=0
00032   Tau=300.
00033   t_src=600.
00034   DO_SOURCE=1
00035   fsmax=0       
00036   RAN_SEED=1
00037 
00038   !setup background
00039   CALL SETUPcoeff(coeff)
00040   do k=1,nz
00041      if(kk(k) .lt. nz_ghostglobal-90) then
00042         FSmax1=sqrt(maxval(&
00043              & (    coeff(:,:,k,15,1)**2+coeff(:,:,k,16,1)**2&
00044              &                                +coeff(:,:,k,17,1)**2)&
00045              &                /(4. *3.1415926)*coeff(:,:,k,12,1)&
00046              & +coeff(:,:,k,5,1) &
00047              & ))
00048      endif
00049      FSmax=max(FSmax,FSmax1)
00050   enddo
00051   
00052   !!set timestep
00053   dt1=min(dx,dy, dz)/FSmax/6.5    !6 with diff*2000
00054   call mpi_reduce (dt1, dt, 1, MPI_REAL, MPI_MIN, 0, MPI_COMM_WORLD, ierr )
00055   dt=15./(int(15./dt)+1) 
00056   call mpi_bcast (dt, 1, MPI_REAL, 0, MPI_COMM_WORLD, ierr)
00057   !full cube sample
00058   nsamp=(300+dt/30.)/(2*dt) 
00059   !slice sample
00060   ssamp=(60+dt/30.)/(2*dt)
00061 
00062   !     INITIALIZE FOURIER TRANSFORMS and derivative arrays and data
00063   call initialize(ider,diff,dt)
00064   call initfftw()
00065   call init_forcing(N_kick)
00066   
00067   !next does IC
00068   
00069   CALL input(ol_old,told)
00070 
00071   if (my_rank .eq. 0) write(*,*) "restarting from step:",ol_old, "t=", t  
00072   t=0
00073   !REPOPULATE RANDOM SOURCE PROPERTIES  
00074   write(*,*) "Repop",ol_old, DO_SOURCE,n_src,N_kick
00075   do ol=0,ol_old-1     
00076      if (DO_SOURCE .eq. 1) then
00077         if (mod(ol,N_kick) .eq. 0) then
00078            ll=mod(ol/N_kick,n_src)+1
00079            write(*,*) ll,t
00080            tt(ll)=t+t_src
00081            x_src=floor(rhc_ran(ran_seed)*nx)
00082            y_src=floor(rhc_ran(ran_seed)*ny) 
00083            do ii=1,nx/2+1
00084               do jj=1,ny
00085                  src_move(ii,jj,ll)=cmplx(cos(x_src*(ii-1)/float(nx)*2*pi),&
00086                       &                            sin(x_src*(ii-1)/float(nx)*2*pi))
00087                  src_move(ii,jj,ll)=src_move(ii,jj,ll)*  &
00088                       &             cmplx(cos(y_src*(jj-1)/float(ny)*2*pi),&
00089                       &                   sin(y_src*(jj-1)/float(ny)*2*pi))
00090               enddo
00091            enddo
00092         endif
00093      endif
00094     t=t+2.*dt
00095   enddo
00096   write(*,*) "t",t,told
00097   t=told
00098 
00099   step=1
00100   call derivs(PE,KE)
00101   mark=-1
00102 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
00103 !!! the time integration LAX-Wendroff !!!
00104 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00105   do ol=ol_old,10000000   
00106      fdold=fd
00107      step=1
00108      call MPI_ghost()
00109      call BCs()
00110      dfddt(:,:,:,1)=centref(fd(:,:,:,4))
00111      dfddt(:,:,:,2)=centref(fd(:,:,:,5))
00112      dfddt(:,:,:,3)=centref(fd(:,:,:,6))
00113      call derivs(PE,KE)
00114 
00115      if (mod(ol,nsamp) .eq. 0) then
00116         mark=ol/nsamp
00117         if(my_rank .eq. 0) write(*,*) "saving cube",t,mark
00118         call outputcube(fd,rd,rpress,rrhoe,rdb,t,mark,ol)
00119      endif
00120      if (mod(ol,ssamp) .eq. 0) then
00121         mark=ol/ssamp
00122         if (my_rank .eq. 0) write(*,*) "saving slice",t,mark,ke,pe
00123         call outputslice(fd,rd,rpress,rrhoe,rdb,t,mark)
00124      endif
00125      if(my_rank .eq. 0) write(*,*) t,PE,KE,ol,nsamp
00126      !END AVERAGE
00127      fd(:,:,1:nz-1,:)=(fd(:,:,1:nz-1,:)+fd(:,:,2:nz  ,:))/2.  
00128 
00129      fd=fd+(dfddt)*dt  
00130      ! MPI PASS STUFF HERE
00131      call MPI_ghost()
00132      call BCs() 
00133      step=2
00134      call derivs(PE,KE)
00135      dfddt(:,:,:,1)=(fdold(:,:,:,4))
00136      dfddt(:,:,:,2)=(fdold(:,:,:,5))
00137      dfddt(:,:,:,3)=(fdold(:,:,:,6))
00138 
00139      !Aaron implemented something vaguely similar - I prefer this
00140      if (DO_SOURCE .eq. 1) then
00141         if (mod(ol,N_kick) .eq. 0) then
00142            ll=mod(ol/N_kick,n_src)+1
00143            tt(ll)=t+t_src
00144            x_src=floor(rhc_ran(ran_seed)*nx)
00145            y_src=floor(rhc_ran(ran_seed)*ny) 
00146            do ii=1,nx/2+1
00147               do jj=1,ny
00148                  src_move(ii,jj,ll)=cmplx(cos(x_src*(ii-1)/float(nx)*2*pi),&
00149                       &                   sin(x_src*(ii-1)/float(nx)*2*pi))
00150                  src_move(ii,jj,ll)=src_move(ii,jj,ll)*  &
00151                       &             cmplx(cos(y_src*(jj-1)/float(ny)*2*pi),&
00152                       &                   sin(y_src*(jj-1)/float(ny)*2*pi))
00153               end do
00154            enddo
00155         endif 
00156         src_fourier=0
00157         do ll=1,n_src 
00158            t_corr=exp(-((t-tt(ll))/(tau))**2)
00159            do ii=1,nx/2+1
00160               do jj=1,ny
00161                  src_fourier(ii,jj,:)=src_fourier(ii,jj,:)+&
00162                       &      fsrc_amp(ii,jj,:)*t_corr*src_move(ii,jj,ll)
00163               end do
00164            end do
00165         enddo
00166         !                                                                               
00167         ! now add the forcing function to the vertical velocity              
00168         dfddt(:,:,:,6)= dfddt(:,:,:,6)+src_fourier
00169      endif
00170      !End random forcing 
00171 
00172      fd=fdold+(dfddt)*2.*dt
00173      t=t+2.*dt
00174      call MPI_ghost()
00175      !===========================
00176      do i1=1,6
00177         temp=d2fdz2(fd(:,:,:,i1))
00178         diffdfddt(:,:,:,i1)=(temp)
00179      enddo
00180      fd(:,:,:,:)=fd(:,:,:,:)+diffdfddt(:,:,:,:)*2*dt*nu*2000 
00181      do i1=1,6
00182         diffdfddt(:,:,:,i1)=fd(:,:,:,i1)*diff(:,:,:)
00183      enddo
00184      fd=fd+diffdfddt*2*dt
00185      call mpi_barrier(MPI_COMM_WORLD,ierr)
00186   enddo
00187   !     =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
00188 end subroutine wave
00189 
00190 !================================================================
00191 
00192 subroutine BCs()
00193   include 'comments.f90'
00194   implicit none
00195   include 'param.f90'
00196   include 'terms3d.f90'
00197   include 'functionlist.f90'
00198   complex  tempin(nx/2+1,ny,2,nvar)
00199   complex tempout(nx/2+1,ny,2,nvar)
00200   integer i     
00201   if(isbottom) then
00202      do i=1,6 
00203         fd(:,:,1,i)=fd(:,:,2,i)*0.99985
00204     enddo
00205   endif   
00206   if(istop) then
00207      fd(:,:,nz,1:6)=(2.*fd(:,:,nz-1,1:6)-fd(:,:,nz-2,1:6))
00208      fd(:,:,nz-1,1:6)=(2.*fd(:,:,nz-1,1:6)-fd(:,:,nz-2,1:6))
00209      !fd(:,:,nz,1:6)=0
00210   endif
00211 end subroutine BCs
00212 !============================================
00213 subroutine MPI_ghost()
00214   include 'comments.f90'
00215   implicit none
00216   include 'param.f90'
00217   include 'terms3d.f90'
00218   include 'functionlist.f90'
00219   complex  tempin(nx/2+1,ny,2,nvar)
00220   complex tempout(nx/2+1,ny,2,nvar)
00221   integer  nmess,up,down,ms
00222   integer*8 i
00223   ms=(nx/2+1)*ny*2*nvar
00224   up=my_rank+1
00225   down=my_rank-1
00226   nmess=12
00227   tempout=fd(:,:,nz-3:nz-2,:)
00228   if(istop .eqv. .false.) then
00229      call mpi_send(tempout,ms,MPI_COMPLEX,&
00230           &                         up,nmess,MPI_COMM_WORLD,ierr)
00231   endif
00232   if(isbottom .eqv. .false.) then
00233              call mpi_recv(tempin ,ms,MPI_COMPLEX,&
00234           &                         down,nmess,MPI_COMM_WORLD,istatus,ierr)
00235              fd(:,:,1:2,:)=tempin
00236   endif 
00237 
00238 
00239   tempout=fd(:,:,3:4,:)
00240   if(isbottom .eqv. .false.) call mpi_send(tempout,ms,MPI_COMPLEX,&
00241           &                         down,nmess,MPI_COMM_WORLD,ierr)
00242   if(istop .eqv. .false.)   then 
00243      call mpi_recv(tempin ,ms,MPI_COMPLEX,&
00244           &                         up,nmess,MPI_COMM_WORLD,istatus,ierr)
00245      fd(:,:,nz-1:nz,:)=tempin
00246   endif  
00247 end subroutine MPI_ghost
00248 
00249 
00250 FUNCTION RHC_RAN (RAN_SEED)
00251   !
00252   ! Random numbers GENERATOR BASED ON NR
00253   !
00254   IMPLICIT NONE
00255   INTEGER RAN_SEED,A,M,Q,R,BM,I
00256   REAL div_M,RHC_RAN
00257   PARAMETER(A=16807,M=2147483647,Q=127773,div_M=1./M,R=2836,BM=123459876)
00258   RAN_SEED=IEOR(RAN_SEED,BM)
00259   I=RAN_SEED/Q
00260   RAN_SEED=A*(RAN_SEED-I*Q)-R*I
00261   IF(RAN_SEED .LT. 0) RAN_SEED=RAN_SEED+M
00262   RHC_RAN=div_m*RAN_SEED
00263   RAN_SEED=IEOR(RAN_SEED,BM)
00264   RETURN
00265 END FUNCTION RHC_RAN
00266 
00267 
00268 
00269