| 1 | subroutine add_noise(u_cplx,v_cplx,efactor, |
|---|
| 2 | * nx,ny,nz,num_dims, |
|---|
| 3 | * wnx,wny,wnz,N,NM1, |
|---|
| 4 | * myid,numprocs,locnx,comm) |
|---|
| 5 | |
|---|
| 6 | c routine to add random noise to the horizontal velocity fields |
|---|
| 7 | c linear vortex component only |
|---|
| 8 | c Energy spectrum defined in function Eofk |
|---|
| 9 | c efact = Enoise/Einput controls the level of the noise |
|---|
| 10 | |
|---|
| 11 | c implicit none |
|---|
| 12 | integer nx,ny,nz,num_dims,N,NM1,nyplanes |
|---|
| 13 | integer myid,numprocs,ierr,locnx,i,j,k |
|---|
| 14 | integer comm |
|---|
| 15 | integer iseed,npts |
|---|
| 16 | parameter (npts=2048) ! >=nx |
|---|
| 17 | double precision urv(npts),vrv(npts),wrv(npts) |
|---|
| 18 | complex u_cplx(locnx,ny+1,nz+1,2),v_cplx(locnx,ny+1,nz+1,2) |
|---|
| 19 | complex w_cplx(locnx,ny+1,nz+1,2) |
|---|
| 20 | real wnx(locnx),wny(ny+1),wnz(nz+1),rmax,efactor |
|---|
| 21 | |
|---|
| 22 | real Ein,Etmp,kappa,xx,utmp,vtmp |
|---|
| 23 | complex zeta |
|---|
| 24 | |
|---|
| 25 | external Enoise |
|---|
| 26 | |
|---|
| 27 | #include "mpif.h" |
|---|
| 28 | |
|---|
| 29 | if( ny .eq. 0 ) then |
|---|
| 30 | rmax = 8./9.*sqrt( 2*wnz( nz/2 )**2 ) ! assumes dx~dz |
|---|
| 31 | else |
|---|
| 32 | rmax = 8./9.*sqrt( 2*wny( ny/2 )**2 ) ! assumes dx~dy |
|---|
| 33 | endif |
|---|
| 34 | |
|---|
| 35 | if( num_dims .eq. 2 ) then |
|---|
| 36 | nyplanes=1 |
|---|
| 37 | elseif( num_dims .eq. 3 ) then |
|---|
| 38 | nyplanes=ny |
|---|
| 39 | endif |
|---|
| 40 | |
|---|
| 41 | c random number initialization |
|---|
| 42 | iseed=myid ! different initial seeds for each processor |
|---|
| 43 | |
|---|
| 44 | c determine HKE of input fields (un-normalized) |
|---|
| 45 | Ein = 0.0 |
|---|
| 46 | do k=1,nz |
|---|
| 47 | do j=1,nyplanes |
|---|
| 48 | do i=1,locnx |
|---|
| 49 | Ein = Ein + .5*( cabs( u_cplx(i,j,k,N) )**2 + cabs( v_cplx(i,j,k,N) )**2 ) |
|---|
| 50 | enddo |
|---|
| 51 | enddo |
|---|
| 52 | enddo |
|---|
| 53 | call MPI_ALLREDUCE(Ein,xx,1,MPI_REAL,MPI_SUM,comm,ierr) ! xx <-- global sum |
|---|
| 54 | Ein=xx |
|---|
| 55 | |
|---|
| 56 | |
|---|
| 57 | c Loop through wavenumber domain |
|---|
| 58 | do k=1,nz |
|---|
| 59 | do j=1,nyplanes |
|---|
| 60 | c obtain npts normal, zero mean, variance=1 |
|---|
| 61 | call zufalli(iseed) ! initialization |
|---|
| 62 | call normalen(npts,urv) |
|---|
| 63 | call normalen(npts,vrv) |
|---|
| 64 | call normalen(npts,wrv) |
|---|
| 65 | |
|---|
| 66 | do i=1,locnx |
|---|
| 67 | |
|---|
| 68 | kappa= sqrt( wnx(i)**2 + wny(j)**2 + wnz(k)**2 )/rmax |
|---|
| 69 | if( kappa .eq. 0 .or. kappa .gt. 1 ) kappa=1.e32 |
|---|
| 70 | xx = Enoise(kappa)/(rmax*kappa) ! i.e. E/sqrt(k^2+l^2+m^2) |
|---|
| 71 | urv(i)=urv(i)*xx |
|---|
| 72 | vrv(i)=vrv(i)*xx |
|---|
| 73 | wrv(i)=wrv(i)*xx |
|---|
| 74 | |
|---|
| 75 | utmp = wny(j)*wrv(i) - wnz(k)*vrv(i) |
|---|
| 76 | vtmp = wnz(k)*urv(i) - wnx(i)*wrv(i) |
|---|
| 77 | |
|---|
| 78 | c keep vortex part only |
|---|
| 79 | zeta = (0.,1.)*wnx(i)*vtmp - (0.,1.)*wny(j)*utmp |
|---|
| 80 | |
|---|
| 81 | c put unscaled noise field in u/v_cplx(:,:,:,NM1) for now |
|---|
| 82 | c set noise mode to zero when kx=ky=0 |
|---|
| 83 | xx=(wnx(i)**2 + wny(j)**2) |
|---|
| 84 | if( xx .gt. 0 ) then |
|---|
| 85 | u_cplx(i,j,k,NM1) = wny(j)*zeta/xx |
|---|
| 86 | v_cplx(i,j,k,NM1) =-wnz(k)*zeta/xx |
|---|
| 87 | else |
|---|
| 88 | u_cplx(i,j,k,NM1) = 0. |
|---|
| 89 | v_cplx(i,j,k,NM1) = 0. |
|---|
| 90 | endif |
|---|
| 91 | |
|---|
| 92 | enddo |
|---|
| 93 | iseed=iseed + numprocs ! get new seed for next set of kx locations |
|---|
| 94 | enddo |
|---|
| 95 | enddo |
|---|
| 96 | |
|---|
| 97 | c determine HKE of noise fields (un-normalized) |
|---|
| 98 | Etmp = 0.0 |
|---|
| 99 | do k=1,nz |
|---|
| 100 | do j=1,nyplanes |
|---|
| 101 | do i=1,locnx |
|---|
| 102 | Etmp = Etmp + .5*( cabs( u_cplx(i,j,k,NM1) )**2 + cabs( v_cplx(i,j,k,NM1) )**2 ) |
|---|
| 103 | enddo |
|---|
| 104 | enddo |
|---|
| 105 | enddo |
|---|
| 106 | call MPI_ALLREDUCE(Etmp,xx,1,MPI_REAL,MPI_SUM,comm,ierr) ! xx <-- global sum |
|---|
| 107 | Etmp=xx |
|---|
| 108 | |
|---|
| 109 | c scale the noise fields to the desired level |
|---|
| 110 | c and add to the input fields |
|---|
| 111 | xx = sqrt(efactor*Ein/Etmp) ! scale amplitudes not energy |
|---|
| 112 | Etmp = 0.0 |
|---|
| 113 | do k=1,nz |
|---|
| 114 | do j=1,nyplanes |
|---|
| 115 | do i=1,locnx |
|---|
| 116 | u_cplx(i,j,k,N) = u_cplx(i,j,k,N) + xx*u_cplx(i,j,k,NM1) |
|---|
| 117 | v_cplx(i,j,k,N) = v_cplx(i,j,k,N) + xx*v_cplx(i,j,k,NM1) |
|---|
| 118 | Etmp = Etmp + .5*( cabs( u_cplx(i,j,k,N) )**2 + cabs( v_cplx(i,j,k,N) )**2 ) |
|---|
| 119 | enddo |
|---|
| 120 | enddo |
|---|
| 121 | enddo |
|---|
| 122 | call MPI_ALLREDUCE(Etmp,xx,1,MPI_REAL,MPI_SUM,comm,ierr) ! xx <-- global sum |
|---|
| 123 | Etmp=xx |
|---|
| 124 | if(myid .eq. 0) write(0,*) ' noise field added to transforms of deterministic ICs ' |
|---|
| 125 | if(myid .eq. 0) write(0,*) ' ratio of total energy to deterministic input HKE: ',Etmp/Ein |
|---|
| 126 | if(myid .eq. 0) write(0,*) ' requested ratio (efactor in problem_params.h) :',efactor |
|---|
| 127 | |
|---|
| 128 | return |
|---|
| 129 | end |
|---|
| 130 | |
|---|
| 131 | real function Enoise(k) |
|---|
| 132 | c |
|---|
| 133 | c define the dimensionless energy spectrum as a function |
|---|
| 134 | c of scaled wavenumber k |
|---|
| 135 | c |
|---|
| 136 | real k |
|---|
| 137 | |
|---|
| 138 | if( k .ge. 1 .or. k .eq. 0 ) then ! outside truncation region or k=0 |
|---|
| 139 | Enoise = 0.0 |
|---|
| 140 | else |
|---|
| 141 | Enoise = 1.0 - k |
|---|
| 142 | endif |
|---|
| 143 | |
|---|
| 144 | return |
|---|
| 145 | end |
|---|
