generateTV Subroutine

  subroutine generateTV()

    integer                      :: i,vx1,vx2,vx3,j
    ! generate random T, typically 300s +/- 75s   
    call rh_init()

    allocate(dtimearr(ntimes))  
    allocate(timearray(ntimes)) 

    if (mype==0) then 
      call randomT(timearray)  
      dtimearr = timearray
      do i = 2, ntimes
        dtimearr(i) = timearray(i) - timearray(i-1)
      enddo
    endif   
    call MPI_BARRIER(icomm,ierrmpi)
    if (npe>1) then
      call MPI_BCAST(timearray,ntimes,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
      call MPI_BCAST(dtimearr,ntimes,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    endif
    

    ! spatial distribution
  
    vx1=ntimes*vlim1
    vx2=ntimes*vlim2
    vx3=ntimes*vlim3
    allocate(va1(vx1))
    allocate(va2(vx2))
    allocate(va3(vx3))
    if (mype==0) then 
       call randomV(1,va1)
       call randomV(2,va2)
       call randomV(3,va3)
    endif
    call MPI_BARRIER(icomm,ierrmpi)
    if (npe>1) then
      call MPI_BCAST(va1,vx1,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
      call MPI_BCAST(va2,vx2,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
      call MPI_BCAST(va3,vx3,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    endif
    

  contains

    subroutine randomT(tb)
        double precision, intent(inout) :: tb(:)

        character(len=100)  :: filename
        double precision                :: tt1,tt,mm
        integer                         :: i

        tt = 0.d0
        do i = 1, ntimes
        call random_number(mm)  
        tt1 = periods + variation * 2. * (mm - 0.5)
        tt = tt + tt1
        tb(i) = tt / unit_time
        end do
    
        write(filename,"(a)") "randomtimes.dat"
        open(unit=21,file=filename,form='formatted',status='replace')
        write(21,'(es12.4)') tb
        close(21)
        
    end subroutine randomT

    subroutine randomV(ival,va)
        integer, intent(in)             :: ival
        double precision, intent(inout) :: va(:)

        character(len=100)  :: filename
        character(len=10)   :: xdirection

        double precision, allocatable   :: vn(:,:),vm(:)
        double precision                :: lambda,ampl,lambda_min,lambda_max,rms,phase
        integer                         :: i, j, kk, vlim
    
        select case (ival)
        case (1)
            vlim = vlim1
        case (2)
            vlim = vlim2
        case (3)
            vlim = vlim3
        end select

        allocate(vn(nwaves, vlim)) 
        allocate(vm(vlim))
        lambda_min = 1.d0/dble(vlim)
        lambda_max = 1.d0
        do kk = 1, ntimes
            do i = 1, nwaves
            lambda = (dble(i-1)/dble(nwaves-1))*(lambda_max-lambda_min)+lambda_min
            ampl = lambda**si 
            call random_number(phase)        
            ! phase to vary between -2pi-->+2pi
            phase=4.0d0*dpi*(phase-0.5d0)
            do j = 1, vlim  
                vn(i,j) = ampl*dsin(2.d0*dpi*(dble(j-1)/dble(vlim-1))/lambda+phase)
            end do
            end do     
            rms = 0.d0
            vm = 0.d0    
            do j = 1, vlim  
            ! sum all the waves in grid point j
            do i = 1, nwaves
                vm(j) = vm(j) + vn(i, j) 
            end do
            rms = rms + vm(j)**2
            end do         
            vm = vm / dsqrt(rms/vlim)    ! normalization        
            do j = 1, vlim
            va(vlim*(kk-1)+j) = vm(j)**2
            end do       
        end do

        write(xdirection,'(I2.2)') ival
        write(filename,"(a)") "randompositions"//TRIM(xdirection)//".dat"
        open(unit=22,file=filename,form='formatted',status='replace')
        write(22,'(es12.4)') va
        close(22)

        deallocate(vm)
        deallocate(vn)

    end subroutine randomV
  end subroutine generateTV