apod Subroutine

   subroutine apod(tab_in, tab_out, long, larg, type_apo, param)
   implicit none
   integer(kind=I4), intent(in )                            :: long     !! *2D array length*
   integer(kind=I4), intent(in )                            :: larg     !! *2D array width*
   character(len=*), intent(in )                            :: type_apo !! *apodization type*
   real   (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab_in   !! *input array*
   real   (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_out  !! *apodized array*
   real   (kind=R8), intent(in ), optional                  :: param    !! *apodized array*

      real   (kind=R8) :: a0, a1, a2, ro, u, v, fct, pun, mun, eps, mi, mj, li, lj, ri, rj, fcti
      integer(kind=I4) :: i, j, dis, ii, jj, lo2, la2

      select case ( type_apo(1:6) )

         case('blackm')
            a0 = 21./50.
            a1 = 01./02.
            a2 = 02./25.
            ri = real(long)  ; rj = real(larg)
            mi = (ri + UN)/2 ; mj = (rj + UN)/2
            li = (ri - UN)/2 ; lj = (rj - UN)/2
            do i = 1, long
            u = (i - mi)/li
               do j = 1, larg
                  v = (j - mj)/lj
                  fct = ( a0 + a1 * cos(PI_R8*u) + a2 * cos(2*PI_R8*u) ) *    &  !
                        ( a0 + a1 * cos(PI_R8*v) + a2 * cos(2*PI_R8*v) )         !
                  tab_out(i, j) = tab_in(i, j)*fct
               enddo
            enddo

         case('tuckey') !http://en.wikipedia.org/wiki/Window_function#Tukey_window
            mun = -UN
            pun = +UN
            eps = 0.25_R8

            if ( present(param) ) eps = param

            dis = nint(eps*long/2)
            tab_out(1:long, 1:larg) = UN

            do i = 0, dis
               ro = 2. * i / ( eps*(long-1) ) - UN
               tab_out(i+1, 1:larg) = 0.5_R8*( UN + cos(PI_R8*ro) )
            enddo
            do i = long - 1 - dis, long - 1
               ro = 2.*i/( eps*(long-1) ) + UN - 2./eps
               tab_out(i+1, 1:larg) = 0.5_R8*( UN + cos(PI_R8*ro) )
            enddo

            dis = nint(eps*larg/2)
            do j = 0, dis
               ro = 2. * j / ( eps*(larg-1) ) - UN
               tab_out(1:long, j+1) = tab_out(1:long, j+1) * 0.5_R8 * ( UN + cos(PI_R8*ro) )
            enddo
            do j = larg-1 -dis, larg-1
               ro = 2. * j / ( eps*(larg-1) ) + UN - 2. / eps
               tab_out(1:long, j+1) = tab_out(1:long, j+1) * 0.5_R8 * ( UN + cos(PI_R8*ro) )
            enddo
            tab_out(1:long, 1:larg) = tab_out(1:long, 1:larg) * tab_in(1:long, 1:larg)

         case('hann__')
            do i = 1, long
               fcti = 0.5 * (1.0 - cos(2.0 * PI_R8 * (i - 1) / (long - 1)))
               do j = 1, larg
                  fct = fcti * 0.5 * (1.0 - cos(2.0 * PI_R8 * (j - 1) / (larg - 1)))
                  tab_out(i, j) = tab_in(i, j)*fct
               enddo
            enddo

         case('welch_')
            ri  = long / 2.000001_R8
            rj  = larg / 2.000001_R8
            lo2 = ceiling( ri )
            la2 = ceiling( rj )
            do ii = lo2 - long, lo2 - 1
               u = (ii / ri)**2
               i = max(ii + lo2, 1)
               do jj = la2 - larg, la2 - 1
                  v = (jj / rj)**2
                  j = max(jj + la2, 1)
                  if ( u + v  > UN ) then
                     tab_out(i, j) = 0
                     cycle
                  endif
                  tab_out(i, j) = tab_in(i, j) * ( 1._R8 - (u + v) )
               enddo
            enddo

         case('no_apo')
            tab_out(1:long, 1:larg) = tab_in(1:long, 1:larg)

         case default
            stop 'apod, apodization type bad choice'

      endselect

   return
   endsubroutine apod