       SUBROUTINE interp_vertical(nlev_init,z_init,p_init
     :         ,nlev_tend,z_tend,p_tend,du_tend,dv_tend,dw_tend
     :         ,dth_tend,drv_tend,dt_tend,dq_tend,
     :         ,z_final,p_final,du_final,dv_final,dw_final
     :         ,dth_final,drv_final,dt_final,dq_final)

       implicit none

c-----------------------------------------------------------------------------------
c Vertical interpolation of tendencies forcing profiles onto initial profiles levels
c-----------------------------------------------------------------------------------
c profil_init.txt nb of levels=nlev_init
c contents:
c z_init, p_init, theta_init, rv_init,u_init,v_init,temp_init,qv_init
c-----------------------------------------------------------------------------------
c profil_tend.txt nb of levels=nlev_tend
c                 nb of forcing times= nbt
c contents:
c z_tend, p_tend, du_tend, dv_tend, dw_tend, dth_tend, drv_tend, dt_tend,dq_tend
c-----------------------------------------------------------------------------------

       integer nlevmax
       parameter (nlevmax=41)
       integer nlev_init,nlev_tend,mxcalc

       real z_init(nlev_init), p_init(nlev_init)
       real z_tend(nlev_tend), p_tend(nlev_tend)
       real du_tend(nlev_tend),dv_tend(nlev_tend)
       real dw_tend(nlev_tend),dth_tend(nlev_tend)
       real drv_tend(nlev_tend),dt_tend(nlev_tend)
       real dq_tend(nlev_tend)

       real z_final(nlev_tend), p_final(nlev_tend)
       real du_final(nlev_tend),dv_final(nlev_tend)
       real dw_final(nlev_tend),dth_final(nlev_tend)
       real drv_final(nlev_tend),dt_final(nlev_tend)
       real dq_final(nlev_tend)

       integer l,k,k1,k2,kp
       real aa,frac,frac1,frac2,fact

       do l = 1, nlev_tend

        if (p_init(l).ge.p_tend(nlev_init)) then

        mxcalc=l
         k1=0
         k2=0

         if (p_init(l).le.p_tend(1)) then

         do k = 1, nlev_init-1
          if (p_init(l).le.p_tend(k)
     :       .and. p_init(l).gt.p_tend(k+1)) then
            k1=k
            k2=k+1
          endif
         enddo

         if (k1.eq.0 .or. k2.eq.0) then
          write(*,*) 'PB! k1, k2 = ',k1,k2
          write(*,*) 'l,p_init(l) = ',l,p_init(l)/100
         do k = 1, nlev_init-1
          write(*,*) 'k,p_tend(k) = ',k,p_tend(k)/100
         enddo
         endif

         frac = (p_tend(k2)-p_init(l))/(p_tend(k2)-p_tend(k1))
         z_final(l)= z_tend(k2) - frac*(z_tend(k2)-z_tend(k1))
         du_final(l)= du_tend(k2) - frac*(du_tend(k2)-du_tend(k1))
         dv_final(l)= dv_tend(k2) - frac*(dv_tend(k2)-dv_tend(k1))
         dw_final(l)= dw_tend(k2) - frac*(dw_tend(k2)-dw_tend(k1))
         dth_final(l)= dth_tend(k2) - frac*(dth_tend(k2)-dth_tend(k1))
         drv_final(l)= drv_tend(k2) - frac*(drv_tend(k2)-drv_tend(k1))
         dt_final(l)= dt_tend(k2) - frac*(dt_tend(k2)-dt_tend(k1))
         dq_final(l)= dq_tend(k2) - frac*(dq_tend(k2)-dq_tend(k1))

         else !p_init>p_tend(1)

         k1=1
         k2=2
         frac1 = (p_init(l)-p_tend(k2))/(p_tend(k1)-p_tend(k2))
         frac2 = (p_init(l)-p_tend(k1))/(p_tend(k1)-p_tend(k2))
         z_final(l)= frac1*z_tend(k1) - frac2*z_tend(k2)
         du_final(l)= frac1*du_tend(k1) - frac2*du_tend(k2)
         dv_final(l)= frac1*dv_tend(k1) - frac2*dv_tend(k2)
         dw_final(l)= frac1*dw_tend(k1) - frac2*dw_tend(k2)
         dth_final(l)= frac1*dth_tend(k1) - frac2*dth_tend(k2)
         dt_final(l)= frac1*dt_tend(k1) - frac2*dt_tend(k2)
         dq_final(l)= frac1*dq_tend(k1) - frac2*dq_tend(k2)

         endif ! p_init.le.p_tend(1)
       endif

       enddo 

       do l = 1,nlev_tend
!      print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
!    $        l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
       enddo

          return
          end
