phi_func Subroutine

    subroutine phi_func(z,L,phim,phih,phim_i,phih_i)

        implicit none

        real, intent(in) :: z,L
        real, intent(out) :: phim,phih,phim_i,phih_i
        real a,b,p,q,pi
        parameter (a=16.,b=5.,p=-0.25,q=-0.5,pi=3.141592653589793)
        real eps

        eps=z/L

        if (eps.ge.0) then
            phim=1.+b*eps
            phim_i=-b*eps
            phih=phim
            phih_i=phim_i
        else
            phim=exp(p*log((1.-a*eps)))
            phih=exp(q*log((1.-a*eps)))
            phim_i=2.*log((1.+1./phim)/2.)+log((1.+1./(phim*phim))/2.)-2.*atan(1./phim)+pi/2.
            phih_i=2.*log((1.+1./phih)/2.)
        endif

    end subroutine phi_func