flow.f90

program flow
use module_globals
use module_realistic
!$ use omp_lib
implicit none
integer(I4B)           :: gsx
integer(I4B)           :: gsy
real(DP),parameter :: pix = 1.0_DP
real(DP),dimension(:,:),allocatable :: surface_dem
integer(I4B),parameter :: infile = 14
integer(I4B),parameter :: outfile = 15
character(len=*),parameter :: outfilename = "flow.dat"
character(len=*),parameter :: infilename = "flow.in"
integer(I4B) :: i,j
real(DP) :: x,y,dx
real(DP) :: dtmp
integer(kind=8) :: recsize
integer(I4B)       :: num_octaves  = 12 ! Number of Perlin noise octaves
real(DP)           :: noise_height = 7e0_DP  ! The magnitude of the noise at the first octave
real(DP)           :: pers = 0.50_DP  ! The relative size scaling at each octave level
real(DP)           :: freq = 2.00_DP     ! Spatial size scale factor multiplier at each octave level
real(DP)           :: xyscale = 0.0100_DP
real(DP),dimension(:,:),allocatable :: anchor,xdistanchor,ydistanchor
real(DP)           :: pdistort_scale = 10._DP
real(DP)           :: xdistort = 0.0_DP,ydistort = 0.0_DP
real(DP)           :: dist_pers = 0.5,dist_freq = 2.0
integer(I4B)       :: distoct = 4
integer(I4B)       :: nseed,frame
integer(I4B), allocatable :: seed(:)
!$ real(DP)        :: t1,t2

! Topographic noise parameters
real(DP) :: xynoise, znoise
real(DP) :: noise,xoff,yoff

!Jordan Turbulence parameters
real(DP)           :: gain1 = 70.0_DP
real(DP)           :: gain  = 0.5_DP
real(DP)           :: warp0 = 0.4_DP
real(DP)           :: warp  = 0.35_DP
real(DP)           :: damp0 = 1.0_DP
real(DP)           :: damp  = 0.8_DP
real(DP)           :: damp_scale  = 1.0_DP
real(DP)           :: jordscale = 0.001_DP

call random_seed(size = nseed)
allocate(seed(nseed))
seed = 42 * (/ (i - 1, i = 1, nseed) /)
call random_seed(put = seed)


open(unit=infile,file=infilename,status='old')
read(infile,*) gsx,gsy
read(infile,*) xoff,yoff
read(infile,*) frame
close(infile)

allocate(anchor(2,num_octaves))
allocate(xdistanchor(2,distoct))
allocate(ydistanchor(2,distoct))
call random_number(anchor)
call random_number(xdistanchor)
call random_number(ydistanchor)


allocate(surface_dem(gsx,gsy))

anchor = anchor * gsx * xyscale

do i = 1,num_octaves
   anchor(1,i) = anchor(1,i) + frame * (i - 1) * xyscale * 0.1
end do


!$ NTHREADS = OMP_get_max_threads() ! In the *parallel* case
!$ write(*,'(a)')      ' OpenMP parameters:'
!$ write(*,'(a)')      ' ------------------'
!$ write(*,'(a,i3,/)') ' Number of threads  = ', nthreads

!$ t1 = omp_get_wtime()
!$omp parallel do default(private) shared(xoff,yoff,xyscale,surface_dem) &
!$omp shared(gsx,gsy,noise_height,freq,pers,num_octaves,anchor) &
!$omp shared(pdistort_scale,dist_freq,dist_pers,distoct,xdistanchor,ydistanchor)&
!$omp shared(gain1,gain,warp0,warp,damp0,damp,damp_scale,jordscale)
do j = 1,gsy
   do i = 1,gsx
      x = (i - 1) * xyscale + xoff
      y = (j - 1) * xyscale + yoff
      xdistort = realistic_turbulence(x, y, xyscale * pdistort_scale, dist_freq, dist_pers, distoct, xdistanchor)
      ydistort = realistic_turbulence(x, y, xyscale * pdistort_scale, dist_freq, dist_pers, distoct, ydistanchor)
      x = x + xdistort
      y = y + ydistort
      noise = 0._DP
      noise = noise + realistic_billowedNoise(x, y, noise_height, freq, pers, num_octaves, anchor)
      noise = noise + realistic_turbulence(x, y, 0.1_DP * noise_height, freq, pers, num_octaves, xyscale + anchor)
      x = x * jordscale / xyscale
      y = y * jordscale / xyscale
      noise = noise + realistic_jordanTurbulence(x,y,freq,gain1,gain,warp0,warp,damp0,damp,damp_scale,num_octaves,&
              2*xyscale + anchor)
      surface_dem(i,j) = noise
   end do
end do
!$omp end parallel do
!$ t2 = omp_get_wtime()
write(*,*) 'minmax: ',minval(surface_dem),maxval(surface_dem)
!$ write(*,*) 'Loop time: ',t2-t1

recsize = sizeof(dtmp) * gsx * gsy
open(outfile,file=outfilename,status='replace',form='unformatted',recl=recsize,access='direct')
write(outfile,rec=1) surface_dem
close(outfile)
deallocate(surface_dem,anchor)
deallocate(xdistanchor,ydistanchor)
end program flow