Subroutine that opens a surface file .sur or .dat
The heights are centred, scaled then put into a vector.
If the scale factor sq is negative, the heights are not scaled when reading .sur
By default, the .sur header is dumped
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| character(len=*), | intent(in) | :: | nom_fic | file name |
||
| real(kind=R8), | intent(in) | :: | mu | if > 0, values are centered |
||
| real(kind=R8), | intent(in) | :: | sq | desired height standard deviation |
||
| real(kind=R8), | intent(out), | dimension(:,:), allocatable | :: | tab_s | height array |
|
| type(SCALE_SURF), | intent(out) | :: | scal | object SCALE_SURF |
subroutine read_surf(nom_fic, mu, sq, tab_s, scal)
implicit none
character(len=*), intent(in) :: nom_fic !! *file name*
real(kind=R8), intent(in) :: mu !! *if > 0, values are centered*
real(kind=R8), intent(in) :: sq !! *desired height standard deviation*
type(SCALE_SURF), intent(out) :: scal !! *object [[SCALE_SURF]]*
real(kind=R8), dimension(:,:), allocatable, intent(out) :: tab_s !! *height array*
integer(kind=I4) :: style, i, ii, j, jj, k, nb, eof, tmp, nx, ny
real(kind=R8) :: sqs, mean, dz, lx, ly, lz
type(OBJ_SURF) :: surf
character(len=3) :: ext
real(kind=R8), dimension(:), allocatable :: x, y, z
i = len_trim(nom_fic)
ext = lower( nom_fic(i-2:i) )
if (ext == 'dat') style = SURF_DAT
if (ext == 'sur') style = SURF_SUR
select case (style)
case (SURF_SUR)
call open_surffile(fichier=trim(nom_fic), surf=surf, scal=scal, dump=.false.)
call trans_surf_tab(surf=surf, tab=tab_s)
case (SURF_DAT)
call get_unit(k)
open(unit=k, file=trim(nom_fic), status='old')
nb = 0
do
read(k, *, iostat=eof)
if ( eof/=0 ) exit
nb = nb +1
enddo
rewind(k)
allocate( x(1:nb) )
allocate( y(1:nb) )
allocate( z(1:nb) )
do i = 1, nb
read(k, *) x(i), y(i), z(i)
enddo
close(k)
! the triplet x, y, z is sorted according x
!-----------------------------------------------------------
call sort_real_2real(g=1, d=nb, rtabref=x(1:nb), rtab1=y(1:nb), rtab2=z(1:nb))
i = 1
do
if ( abs(x(i) -x(1)) > 100*EPS_R8 ) exit
i = i +1
enddo
scal%dx = abs(x(i) -x(1))
j = 1
do
if (abs(x(j +1) -x(j))>1.0e-10) exit
j = j +1
enddo
ny = j ! number of same abscissae for a given column
if (mod(nb, ny) /=0 ) STOP 'READ_SURF, non rectangular mesh'
nx = nb/ny
scal%xres = nx
scal%yres = ny
do i = 1, nx
ii = (i-1)*ny +1
jj = ii +ny -1
call sort_real_1real(g=1, d=ny, rtabref=y(ii:jj), rtab1=z(ii:jj))
enddo
j = 1
do
if ( abs(y(j) -y(1)) > 100*EPS_R8 ) exit
j = j +1
enddo
scal%dy = abs(y(j) -y(1))
allocate( tab_s(1:nx, 1:ny) )
k = 0
do i = 1, nx
do j = 1, ny
k = k +1
tab_s(i, j) = z(k)
enddo
enddo
lx = maxval(x) -minval(x)
ly = maxval(y) -minval(y)
lz = maxval(tab_s) -minval(tab_s)
scal%dz = lz/(nx*ny)
scal%xlength_unit = 'm ' ; scal%dx_unit = 'm '
scal%ylength_unit = 'm ' ; scal%dy_unit = 'm '
scal%zlength_unit = 'm ' ; scal%dz_unit = 'm '
!call sort_real(g=1, d=nx*ny, rtabref=z(1:nx*ny))
!scal%dz = 1.e+10
!do i = 2, nx*ny
! dz = z(i) -z(i-1)
! if (abs(dz) < 100*EPS_R8) cycle
! scal%dz = min(scal%dz, dz)
!enddo
!scal%dz = dz
deallocate(x, y, z)
endselect
nx = scal%xres
ny = scal%yres
! centering ?
if (mu > 0.) then
mean = sum(tab_s(1:nx, 1:ny)) / (nx * ny)
scal%mu = mean
tab_s(1:nx, 1:ny) = tab_s(1:nx, 1:ny) -mean
endif
! scaling ?
if (sq > 0.) then
sqs = (sum(tab_s(1:nx, 1:ny) ** 2) / (nx * ny)) ** (0.5_R8)
scal%si = sqs
tab_s(1:nx, 1:ny) = tab_s(1:nx, 1:ny) * sq / sqs
endif
return
endsubroutine read_surf