program test_morpho
use data_arch, only : I4, R8, PI_R8
use miscellaneous, only : get_unit
use surfile, only : init_scal, read_surf, write_surf, SCALE_SURF
use sort_arrays, only : sort_array2
use stat_mom, only : calc_moments, moment_stat
use morpho, only : calcul_normales, surf_area, count_cell, erode_dilate, def_masque
implicit none
real (kind=R8), allocatable, dimension(:,:) :: mask, array
real (kind=R8), allocatable, dimension(:) :: vec, vec_tmp
integer(kind=I4), allocatable, dimension(:,:) :: imask
type(SCALE_SURF) :: scal_surf
integer(kind=I4) :: i, nnx, nny, nb_cells
real (kind=R8) :: median_size, c1, c2, topo, mintab, maxtab
integer(kind=I4) :: ni, nj, ii, jj, j, k, ua
real (kind=R8) :: angle, dx, dy, t, cp, sp, res_hori, area
type(moment_stat) :: mom
!----------------------------------------------------------------------------------
!---------- Returns the number of cells of a given mask ---------------------------
!----------------------------------------------------------------------------------
call read_surf( nom_fic = "sur/mask.sur", & !
tab_s = mask, & !
scal = scal_surf ) !
nnx = scal_surf%xres
nny = scal_surf%yres
allocate( imask(1:nnx, 1:nny) )
imask(1:nnx, 1:nny) = nint( mask(1:nnx, 1:nny) )
call count_cell( msk = imask(1:nnx, 1:nny), & !
long = nnx, & !
larg = nny, & !
nbr_cell = nb_cells, & !
med_cell = median_size ) !
call write_surf( nom_fic = "out/res_count_cell.sur", & !
tab_s = real( imask(1:nnx, 1:nny), kind = R8 ), & !
scal = scal_surf ) !
write(*, '(a, T25, I3, a, E12.4, a)') 'initial mask:', nb_cells, ' cells, ', median_size, ' % surface'
!----------------------------------------------------------------------------------
!---- Returns the number of cells of a given mask after erosion -------------------
!----------------------------------------------------------------------------------
imask(1:nnx, 1:nny) = nint( mask(1:nnx, 1:nny) )
call erode_dilate( msk = imask(1:nnx, 1:nny), & !
long = nnx, & !
larg = nny, & !
val = 6, & !
act = 'erode' ) !
call count_cell( msk = imask(1:nnx, 1:nny), & !
long = nnx, & !
larg = nny, & !
nbr_cell = nb_cells, & !
med_cell = median_size ) !
call write_surf( nom_fic = "out/res_count_cell_erode.sur", & !
tab_s = real( imask(1:nnx, 1:nny), kind = R8 ), & !
scal = scal_surf ) !
write(*, '(a, T25, I3, a, E12.4, a)') 'mask after erode:', nb_cells, ' cells, ', median_size, ' % surface'
!----------------------------------------------------------------------------------
!----- Returns the number of cells of a given mask after erosion and dilation -----
!----------------------------------------------------------------------------------
where ( imask(1:nnx, 1:nny) >= 1 ) imask(1:nnx, 1:nny) = 1
call erode_dilate( msk = imask(1:nnx, 1:nny), & !
long = nnx, & !
larg = nny, & !
val = 6, & !
act = 'dilat' ) !
call count_cell( msk = imask(1:nnx, 1:nny), & !
long = nnx, & !
larg = nny, & !
nbr_cell = nb_cells, & !
med_cell = median_size ) !
call write_surf( nom_fic = "out/res_count_cell_erode_dilat.sur", & !
tab_s = real( imask(1:nnx, 1:nny), kind = R8 ), & !
scal = scal_surf ) !
write(*, '(a, T25, I3, a, E12.4, a)') 'mask erode + dilate:', nb_cells, ' cells, ', median_size, ' % surface'
write(*,*) '--------------------------------------------------------'
!----------------------------------------------------------------------------------
!---------- Returns the % age of heights above c2, once c1 heights are removed ----
!----------------------------------------------------------------------------------
allocate( array( 1:nnx, 1:nny ) )
allocate( vec( 1:nnx * nny ), vec_tmp( 1:nnx * nny ) )
! thresholds
c1 = 0.15_R8
c2 = 0.85_R8
! surface made of increasing integers
vec( 1:nnx * nny ) = [ ( i, i = 1, nnx * nny ) ]
mintab = minval( vec(1:nnx*nny) )
maxtab = maxval( vec(1:nnx*nny) )
! heights rescaled between 0 and 1
vec_tmp(1:nnx*nny) = ( vec(1:nnx*nny) - mintab )/( maxtab - mintab )
! shuffle heights
call random_number( vec(1:nnx*nny) )
call sort_array2(tab_inout = vec(1:nnx*nny), & !
tab1 = vec_tmp(1:nnx*nny), & !
n = nnx*nny) !
array( 1:nnx, 1:nny ) = reshape( vec_tmp(1:nnx*nny), [nnx, nny] )
! mask heights above c2, when heights below c1 are ignored
call def_masque( msk = imask(1:nnx, 1:nny), & !
tab = array(1:nnx, 1:nny), & !
long = nnx , & !
larg = nny , & !
crit1 = c1, & !
crit2 = c2, & !
top = topo ) !
call write_surf( nom_fic = "out/mask_after_threshold.sur", & !
tab_s = real( imask(1:nnx, 1:nny), kind = R8 ), & !
scal = scal_surf ) !
! topo is the fraction of heights that are involved
write(*,*) 'percentage of heights abobe 85%, without the first 15%: ', topo ! = 100. - 100 * ( c1 + c2 * (1. - c1) )
write(*,*) '--------------------------------------------------------'
deallocate( array, mask, imask, vec, vec_tmp )
!----------------------------------------------------------------------------------
!--------- Returns the % age of surface nearly horizontal -------------------------
!----------------------------------------------------------------------------------
nnx = 512 ; nny = 512
allocate( array( 1:nnx, 1:nny ) )
allocate( vec( 1:nnx * nny ) )
! the domain is divided into 4x4 patches
ni = 512 / 4
nj = 512 / 4
! cone half angle
angle = 5._R8
! lag along x and y
dx = 1.e-7_R8
dy = 1.e-7_R8
k = 0
do i = 1, 4
do j = 1, 4
k = k + 1 ! subdomain number
! the subdomain is a plane of equation: z = -a.x -b.y -c
! the normal coordinates are (-a, -b, +1)/sqrt(a**2 + b**2 + 1)
! and the scalar product with \vec{k} is 1/sqrt(a**2 + b**2 + 1)
! which is also cos(\theta).
! Therefore a**2 + b**2 = tan(\theta)**2, so if a = cos(\phi).tan(\theta)
! b = sin(\phi).tan(\theta), the relationship is satisfied (whatever \phi).
! With \theta = k * 0.11 * angle, 9 subdomains are concerned.
t = tan( k * 0.11 * angle * PI_R8 / 180 )
call random_number( cp )
cp = 2 * (cp - 0.5)
call random_number( sp )
if (sp > 0.5_R8) then
sp = + sqrt( 1._R8 - cp**2)
else
sp = - sqrt( 1._R8 - cp**2)
endif
do ii = (i - 1) * ni + 1, i * ni
do jj = (j - 1) * nj + 1, j * nj
array(ii, jj) = plane( a = cp * t, & !
b = sp * t, & !
c = 0._R8, & !
x = ii * dx, & !
y = jj * dy ) !
enddo
enddo
enddo
enddo
call calc_moments( tab = reshape( array(1:nnx, 1:nny), [nnx * nny] ), & !
mx = mom, & !
nb_mom = 2 ) !
array(1:nnx, 1:nny) = ( array(1:nnx, 1:nny) - mom%mu ) / mom%si
call init_scal( scal = scal_surf, & !
nx = nnx, & !
ny = nny, & !
lx = nnx * dx, & !
ly = nny * dy, & !
unit_z = 'm ' ) !
call write_surf( nom_fic = "out/hori.sur", & !
tab_s = array(1:nnx, 1:nny), & !
scal = scal_surf ) !
call calcul_normales( tab_in = array(1:nnx, 1:nny), & !
long = nnx, & !
larg = nny, & !
scale_xyz = [dx, dy, mom%si], & !
cone_angle = angle, & !
print_mask = .true., & !
hori = res_hori ) !
call surf_area( tab_in = array(1:nnx, 1:nny), & !
long = nnx, & !
larg = nny, & !
scale_xyz = [dx, dy, mom%si], & !
aire = area ) !
call get_unit(ua)
open( unit = ua, file = "out/mask_angle.txt")
read( ua, * ) ( vec(i), i = 1, nnx * nny )
close( ua )
call write_surf( nom_fic = "out/mask_hori.sur", & !
tab_s = reshape( vec(1:nnx*nny), [nnx, nny] ), & !
scal = scal_surf ) !
write(*, *) 'Percentage of nearly horizontal surface (+/- 5°): ', res_hori
write(*, *) 'Relativea area minus 1 and z scale: ', area, mom%si
deallocate( array, vec )
contains
real(kind=R8) function plane(a, b, c, x, y)
implicit none
real(kind=R8), intent(in) :: a, b, c, x, y
plane = a * x + b * y + c
return
endfunction plane
endprogram test_morpho