High performance exponential of a skew Hermitian matrix in Fortran 95
Your code looks good. When s is small, there should be no need to use BLAS instead of matmul(). I like you exploit the shorthand array notation and the forall construct.

Debatable is your use of explicit shaped arrays. They could be very slightly faster, because they are contiguous, but in many cases the compiler will have to copy the arrays when calling your routines and the overall code will be slower. The preferred way in modern code are assumed shape arrays, e.g.

 complex(complex_kind), intent(in) :: H(:,:)
 complex(complex_kind)             :: B(1:ubound(H,1),1:ubound(H,1))


but there may be reasons to use your way and it may be shorter.

Generally avoid using star notation real*4 *8 *16 for sizes of your variables. It is non-standard and obsolete. Use real(some_kind_constant).

You can use selected_real_kind() and selected_int_kind() (now preferred), to get the constants.

Or if you don't mind Fortran 2008 features and need to know the size in bits you can use kind constants from the iso_fortran_env module, like integer(int32), real(real64) and so on. 

Another possibility from Fortran 2003 is to use kind constants from the iso_c_binding module, like integer(c_int) if you need interoperability with C code.