read_surf Subroutine

   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 ), optional  :: mu                    !! *desired mean*
   real(kind=R8),    intent(in ), optional  :: 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_array2(tab_inout = x(1:nb),           &  !
                                  tab1 = y(1:nb),           &  !
                                  tab2 = z(1:nb), n = 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_array2(tab_inout = y(ii:jj),                    &  !
                                     tab1 = z(ii:jj), n = jj - ii + 1)      !
            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

      scal%mu =  sum(  tab_s(1:nx, 1:ny)               ) / (nx * ny)
      scal%si = (sum( (tab_s(1:nx, 1:ny) -scal%mu) ** 2) / (nx * ny)) ** (0.5_R8)

      ! centering ?
      if ( present(mu) ) then
         tab_s(1:nx, 1:ny) = (tab_s(1:nx, 1:ny) -scal%mu) + mu
      endif

      ! scaling ?
      if ( present(sq) ) then
         mean =  sum(  tab_s(1:nx, 1:ny)            ) / (nx * ny)
         sqs  = (sum( (tab_s(1:nx, 1:ny) -mean) ** 2) / (nx * ny)) ** (0.5_R8)

         tab_s(1:nx, 1:ny) = sq*(tab_s(1:nx, 1:ny) -mean)/sqs + mean
      endif

   return
   endsubroutine read_surf