C CONJP2    SOURCE    CB215821  16/04/21    21:15:57     8920
       SUBROUTINE CONJP2(JL,JR,WVEC_L,WVEC_R,NVECT,TVECT,
     &                   ga)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  CONJP2
C
C DESCRIPTION       :  Voir KONJP2
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  S. KUDRIAKOV, DM2S/SFME/LTMF
C
C************************************************************************
C
c----------------------------------------------------------------------
c GENERAL DESCRIPTION:
c   This subroutine provides the jacobians which are the derivatives
c   of the numerical flux function defined at the cell interface
c   with respect to the primitive variables of the left and right
c   cells (relative to the cell interface).
c
c EQUATIONS: 2D Euler equations of gas dynamics
c
c
c REFERENCE: JCP, 129, 364-382 (1996)
c            " A Sequel to AUSM: AUSM+ ".
c----------------------------------------------------------------------
c INPUT:
c
c      alpha   -- parameter of the AUSM+ scheme in the Pressure function;
c                 ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
c
c      beta    -- parameter of the AUSM+ scheme in the Mach function;
c                 ( -1/16 <= beta <= 1/2 )  (imposed as a parameter)
c
c      wvec_l  -- vector of the primitive variables (rho,ux,uy,p) at the
c                 left cell;
c
c      wvec_r  -- vector of the primitive variables (rho,ux,uy,p) at the
c                 right cell;
c
c      nvect   -- normal vector to the interface (2 components in 2D);
c
c      tvect   -- tangential vector to the interface;
c
c      ga      -- ratio of the specific heats (assumed constant)
c----------------------------------------------------------------------
c
c OUTPUT:
c
c      jl     -- jakobian matrix 4 by 4 - derivatives of the numerical
c                 flux function with respect to the primitive variables
c                 from the left cell;
c
c      jr     -- jakobian matrix 4 by 4 - derivatives of the numerical
c                 flux function with respect to the primitive variables
c                 from the right cell.
c
c N.B. Numerical flux function is given by
c      \rho un
c      \rho un ux + p (n.x)
c      \rho un uy + p (n.y)
c      \rho un ht
c----------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
       real*8 wvec_l(4),wvec_r(4)
       real*8 nvect(2),tvect(2)
       real*8 jl(4,4),jr(4,4),f(4)
       real*8 alpha,beta
       real*8 ga, gm1,temph
       real*8 n_x,n_y,t_x,t_y
       real*8 un_l, un_r, ut_l, ut_r
       real*8 ml,mr,Mplus,Mmin,mmid
       real*8 mpl_m, mmin_m,am
       real*8 rold_l,uold_l,vold_l,pold_l,eold_l
       real*8 rold_r,uold_r,vold_r,pold_r,eold_r
       real*8 Pplus,Pmin,pmid
       real*8 hr_l,hr_r,det
       real*8 br1,br2,br3,br4,temp_l,temp_r,brac_l,brac_r
       real*8 aleft, arigh
       real*8 damr_l,damr_r,damu_l,damu_r
       real*8 damv_l,damv_r,damp_l,damp_r
       real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
       real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
       real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
       real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
       real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
       real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
       real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
       real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
       real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
       real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
       real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
       real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
       real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
       real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
       real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
       real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
       real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
       real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
       real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
       real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
       parameter(alpha = 0.1875D0, beta = 0.125D0)
c-----------------------
       gm1=ga-1.0d0
c-----------------------
       n_x=nvect(1)
       n_y=nvect(2)
       t_x=tvect(1)
       t_y=tvect(2)
c----------------------------
c----------------------------
       rold_l=wvec_l(1)
       uold_l=wvec_l(2)
       vold_l=wvec_l(3)
       pold_l=wvec_l(4)
c-----------------------
       rold_r=wvec_r(1)
       uold_r=wvec_r(2)
       vold_r=wvec_r(3)
       pold_r=wvec_r(4)
c------------------------------------------------------------------
c Computation of the specific total energy on the left and right.
c------------------------------------------------------------------
       eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
       eold_l=eold_l+pold_l/(gm1*rold_l)
       eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
       eold_r=eold_r+pold_r/(gm1*rold_r)
c-------------------------------------------------------------------
c Computation of the speed of sound and its derivatives;
c numerical speed of sound at the interface is taken as an average
c of the speeds of sounds of the neighbouring cells
c-------------------------------------------------------------------
       aleft=SQRT(ga*pold_l/rold_l)
       arigh=SQRT(ga*pold_r/rold_r)
       am=0.5d0*(aleft+arigh)
c-------------------------------------------------------------------
         damr_r=-arigh/(4.0d0*rold_r)
         damu_r=0.0d0
         damv_r=0.0d0
         damp_r=ga/(4.0d0*arigh*rold_r)
c-----------------------
         damr_l=-aleft/(4.0d0*rold_l)
         damu_l=0.0d0
         damv_l=0.0d0
         damp_l=ga/(4.0d0*aleft*rold_l)
c-------------------------------------------------------------------
c Computing numerical Mach number and its derivatives,
c see p.370, under (A1).
c-------------------------------------------------------------------
       un_l=uold_l*n_x+vold_l*n_y
       un_r=uold_r*n_x+vold_r*n_y
       ut_l=uold_l*t_x+vold_l*t_y
       ut_r=uold_r*t_x+vold_r*t_y
c-------------------------------------------------------------------
       ml=un_l/am
       mr=un_r/am
c-------------------------------------------------------------------
c  Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
c  see (19a) and (19b), p.367.
c-------------------------------------------------------------------
       if(ABS(ml) .ge. 1.0d0) then
         Mplus=(ml+ABS(ml))/2.0d0
       else
         Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
         Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
       endif
c-------------------------------------------------------------------
       if(ABS(mr) .ge. 1.0d0) then
         Mmin=(mr-ABS(mr))/2.0d0
       else
         Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
         Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
       endif
c-------------------------------------------------------------------
c Derivatives of ml and mr with respect to both sets of primitive
c variables: left  and right.
c-------------------------------------------------------------------
       temp_l=-un_l/(am*am)
       temp_r=-un_r/(am*am)
c--------
       dmlr_l=temp_l*damr_l
       dmlr_r=temp_l*damr_r
       dmrr_l=temp_r*damr_l
       dmrr_r=temp_r*damr_r
c--------
       dmlu_l=n_x/am+temp_l*damu_l
       dmlu_r=temp_l*damu_r
       dmru_l=temp_r*damu_l
       dmru_r=n_x/am+temp_r*damu_r
c--------
       dmlv_l=n_y/am+temp_l*damv_l
       dmlv_r=temp_l*damv_r
       dmrv_l=temp_r*damv_l
       dmrv_r=n_y/am+temp_r*damv_r
c--------
       dmlp_l=temp_l*damp_l
       dmlp_r=temp_l*damp_r
       dmrp_l=temp_r*damp_l
       dmrp_r=temp_r*damp_r
c-----------------------------------------------------------
c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
c-----------------------------------------------------------
       mmid=Mplus+Mmin
c-----------------------------------------------------------
c Computing the derivatives of M+ and M-
c-----------------------------------------------------------
       if(ml .ge. 1.0d0) then
         dMpr_l=dmlr_l
         dMpu_l=dmlu_l
         dMpv_l=dmlv_l
         dMpp_l=dmlp_l
c--------------------
         dMpr_r=dmlr_r
         dMpu_r=dmlu_r
         dMpv_r=dmlv_r
         dMpp_r=dmlp_r
       else
        if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
         temph=(ml+1.0d0)/2.0d0
         dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
         dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
         dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
         dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
c--------------------
           dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
           dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
           dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
           dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
         else
           dMpr_l=0.0d0
           dMpu_l=0.0d0
           dMpv_l=0.0d0
           dMpp_l=0.0d0
c---------------------
           dMpr_r=0.0d0
           dMpu_r=0.0d0
           dMpv_r=0.0d0
           dMpp_r=0.0d0
         endif
       endif
c-----------------------------------------------------------
       if(mr .ge. 1.0d0) then
         dMmr_l=0.0d0
         dMmu_l=0.0d0
         dMmv_l=0.0d0
         dMmp_l=0.0d0
c---------------------
         dMmr_r=0.0d0
         dMmu_r=0.0d0
         dMmv_r=0.0d0
         dMmp_r=0.0d0
       else
         if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
           temph=(-mr+1.0d0)/2.0d0
           dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
           dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
           dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
           dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
c--------------------
           dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
           dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
           dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
           dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
         else
           dMmr_l=dmrr_l
           dMmu_l=dmru_l
           dMmv_l=dmrv_l
           dMmp_l=dmrp_l
c--------------------
           dMmr_r=dmrr_r
           dMmu_r=dmru_r
           dMmv_r=dmrv_r
           dMmp_r=dmrp_r
         endif
       endif
c-----------------------------------------------------------------
c computing the derivatives of m_{1/2} (notation as in the paper)
c-----------------------------------------------------------------
       dmir_l=dMpr_l+dMmr_l
       dmir_r=dMpr_r+dMmr_r
c-------------
       dmiu_l=dMpu_l+dMmu_l
       dmiu_r=dMpu_r+dMmu_r
c-------------
       dmiv_l=dMpv_l+dMmv_l
       dmiv_r=dMpv_r+dMmv_r
c-------------
       dmip_l=dMpp_l+dMmp_l
       dmip_r=dMpp_r+dMmp_r
c----------------------------------------------------------------
c computing the main convective variables and their derivatives
c mpl_m is m^{+}_{1/2} (paper's notation) and
c mmin_m  is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
c----------------------------------------------------------------
       if(mmid .ge. 0.0d0) then
         mpl_m = mmid
         d2mr_l=dmir_l
         d2mu_l=dmiu_l
         d2mv_l=dmiv_l
         d2mp_l=dmip_l
c------------
         d2mr_r=dmir_r
         d2mu_r=dmiu_r
         d2mv_r=dmiv_r
         d2mp_r=dmip_r
c------------
       else
         mpl_m = 0.0d0
         d2mr_l=0.0d0
         d2mu_l=0.0d0
         d2mv_l=0.0d0
         d2mp_l=0.0d0
c------------
         d2mr_r=0.0d0
         d2mu_r=0.0d0
         d2mv_r=0.0d0
         d2mp_r=0.0d0
       endif
c---------------------------------------------
       if(mmid .le. 0.0d0) then
         mmin_m = mmid
         d3mr_l=dmir_l
         d3mu_l=dmiu_l
         d3mv_l=dmiv_l
         d3mp_l=dmip_l
c------------
         d3mr_r=dmir_r
         d3mu_r=dmiu_r
         d3mv_r=dmiv_r
         d3mp_r=dmip_r
c------------
       else
         mmin_m = 0.0d0
         d3mr_l=0.0d0
         d3mu_l=0.0d0
         d3mv_l=0.0d0
         d3mp_l=0.0d0
c------------
         d3mr_r=0.0d0
         d3mu_r=0.0d0
         d3mv_r=0.0d0
         d3mp_r=0.0d0
       endif
c---------------------------------------------------------------
c Computing the calligraphic P+ and P- with their derivatives,
c see (21a) & (21b) on p.368.
c---------------------------------------------------------------
       if(ml .ge. 1.0d0) then
         Pplus = 1.0d0
       else
         if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
           Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
           Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
         else
           Pplus = 0.0d0
         endif
       endif
c---------------------------------------------------------------
       if(mr .ge. 1.0d0) then
         Pmin = 0.0d0
       else
         if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
           Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
           Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
         else
           Pmin = 1.0d0
         endif
       endif
c---------------------------------------------------------------
       brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
       brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
       brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
c--------------
       brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
       brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
       brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
c---------------------------------------------------------------
       if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
         dPpr_l=brac_l*dmlr_l
         dPpr_r=brac_l*dmlr_r
c------------
         dPpu_l=brac_l*dmlu_l
         dPpu_r=brac_l*dmlu_r
c------------
         dPpv_l=brac_l*dmlv_l
         dPpv_r=brac_l*dmlv_r
c------------
         dPpp_l=brac_l*dmlp_l
         dPpp_r=brac_l*dmlp_r
c------------
       else
         dPpr_l=0.0d0
         dPpr_r=0.0d0
c-----------
         dPpu_l=0.0d0
         dPpu_r=0.0d0
c-----------
         dPpv_l=0.0d0
         dPpv_r=0.0d0
c-----------
         dPpp_l=0.0d0
         dPpp_r=0.0d0
c-----------
       endif
c---------------------------------------------------------------
       if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
         dPmr_l=brac_r*dmrr_l
         dPmr_r=brac_r*dmrr_r
c------------
         dPmu_l=brac_r*dmru_l
         dPmu_r=brac_r*dmru_r
c------------
         dPmv_l=brac_r*dmrv_l
         dPmv_r=brac_r*dmrv_r
c------------
         dPmp_l=brac_r*dmrp_l
         dPmp_r=brac_r*dmrp_r
c------------
       else
         dPmr_l=0.0d0
         dPmr_r=0.0d0
c-----------
         dPmu_l=0.0d0
         dPmu_r=0.0d0
c-----------
         dPmv_l=0.0d0
         dPmv_r=0.0d0
c-----------
         dPmp_l=0.0d0
         dPmp_r=0.0d0
c-----------
       endif
c-------------------------------------------------------------------
c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
c-------------------------------------------------------------------
       pmid=Pplus*pold_l + Pmin*pold_r
       dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
       dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
       dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
       dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
c----------------------------
       dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
       dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
       dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
       dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
c---------------------------------------------------------------------
c computing JACOBIAN as a derivative of the numerical flux function
c with respect to the primitive variables.
c Notation: jl(i,j) --- is the derivative of the i-component of the
c           flux function with respect to the j-component of the
c           vector of primitive variables of the left state.
c           jr(i,j) --- is the derivative of the i-component of the
c           flux function with respect to the j-component of the
c           vector of primitive variables of the right state.
c---------------------------------------------------------------------
       f(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
c---------------------------------------------------------------------
       jl(1,1)=damr_l*f(1)/am+am*(d2mr_l*rold_l+mpl_m)
       jl(1,1)=jl(1,1)+am*d3mr_l*rold_r
       jl(1,2)=damu_l*f(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
       jl(1,3)=damv_l*f(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
       jl(1,4)=damp_l*f(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
c------------------------------------
       jr(1,1)=damr_r*f(1)/am+am*(d2mr_r*rold_l+mmin_m)
       jr(1,1)=jr(1,1)+am*d3mr_r*rold_r
       jr(1,2)=damu_r*f(1)/am+am*(d2mu_r*rold_l+d3mu_r*rold_r)
       jr(1,3)=damv_r*f(1)/am+am*(d2mv_r*rold_l+d3mv_r*rold_r)
       jr(1,4)=damp_r*f(1)/am+am*(d2mp_r*rold_l+d3mp_r*rold_r)
c------------------------------------
c------------------------------------
       br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
       br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
       f(2)=n_x*(am*br1+pmid)+am*t_x*br2
c------------------
       det=n_x*t_y-n_y*t_x
c---------------------------------------------------------
       br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
       br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
       jl(2,1)=n_x*(damr_l*br1+am*br3+dpir_l)
       jl(2,1)=jl(2,1)+t_x*damr_l*br2+am*t_x*br4
c-------------------
       br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
       br3=br3+d3mu_l*rold_r*un_r
       br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y)/det
       br4=br4+d3mu_l*rold_r*ut_r
       jl(2,2)=n_x*(damu_l*br1+am*br3+dpiu_l)
       jl(2,2)=jl(2,2)+t_x*damu_l*br2+am*t_x*br4
c-------------------
       br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x)/det
       br3=br3+d3mv_l*rold_r*un_r
       br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
       br4=br4+d3mv_l*rold_r*ut_r
       jl(2,3)=n_x*(damv_l*br1+am*br3+dpiv_l)
       jl(2,3)=jl(2,3)+t_x*damv_l*br2+am*t_x*br4
c-------------------
       br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
       br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
       jl(2,4)=n_x*(damp_l*br1+am*br3+dpip_l)
       jl(2,4)=jl(2,4)+t_x*damp_l*br2+am*t_x*br4
c-------------------------------------------------------------
       br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
       br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
       jr(2,1)=n_x*(damr_r*br1+am*br3+dpir_r)
       jr(2,1)=jr(2,1)+t_x*damr_r*br2+am*t_x*br4
c-------------------
       br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
       br3=br3+d3mu_r*rold_r*un_r
       br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
       br4=br4+d3mu_r*rold_r*ut_r
       jr(2,2)=n_x*(damu_r*br1+am*br3+dpiu_r)
       jr(2,2)=jr(2,2)+t_x*damu_r*br2+am*t_x*br4
c-------------------
       br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
       br3=br3+d3mv_r*rold_r*un_r
       br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
       br4=br4+d3mv_r*rold_r*ut_r
       jr(2,3)=n_x*(damv_r*br1+am*br3+dpiv_r)
       jr(2,3)=jr(2,3)+t_x*damv_r*br2+am*t_x*br4
c-------------------
       br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
       br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
       jr(2,4)=n_x*(damp_r*br1+am*br3+dpip_r)
       jr(2,4)=jr(2,4)+t_x*damp_r*br2+am*t_x*br4
c-------------------------------------------------------------
c-------------------------------------------------------------
       br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
       br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
       f(3)=n_y*(am*br1+pmid)+am*t_y*br2
       br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
       br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
       jl(3,1)=n_y*(damr_l*br1+am*br3+dpir_l)
       jl(3,1)=jl(3,1)+t_y*damr_l*br2+am*t_y*br4
c-------------------
       br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
       br3=br3+d3mu_l*rold_r*un_r
       br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y/det)
       br4=br4+d3mu_l*rold_r*ut_r
       jl(3,2)=n_y*(damu_l*br1+am*br3+dpiu_l)
       jl(3,2)=jl(3,2)+t_y*damu_l*br2+am*t_y*br4
c-------------------
       br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x/det)
       br3=br3+d3mv_l*rold_r*un_r
       br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
       br4=br4+d3mv_l*rold_r*ut_r
       jl(3,3)=n_y*(damv_l*br1+am*br3+dpiv_l)
       jl(3,3)=jl(3,3)+t_y*damv_l*br2+am*t_y*br4
c-------------------
       br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
       br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
       jl(3,4)=n_y*(damp_l*br1+am*br3+dpip_l)
       jl(3,4)=jl(3,4)+t_y*damp_l*br2+am*t_y*br4
c-------------------------------------------------------------
       br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
       br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
       jr(3,1)=n_y*(damr_r*br1+am*br3+dpir_r)
       jr(3,1)=jr(3,1)+t_y*damr_r*br2+am*t_y*br4
c-------------------
       br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
       br3=br3+d3mu_r*rold_r*un_r
       br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
       br4=br4+d3mu_r*rold_r*ut_r
       jr(3,2)=n_y*(damu_r*br1+am*br3+dpiu_r)
       jr(3,2)=jr(3,2)+t_y*damu_r*br2+am*t_y*br4
c-------------------
       br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
       br3=br3+d3mv_r*rold_r*un_r
       br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
       br4=br4+d3mv_r*rold_r*ut_r
       jr(3,3)=n_y*(damv_r*br1+am*br3+dpiv_r)
       jr(3,3)=jr(3,3)+t_y*damv_r*br2+am*t_y*br4
c-------------------
       br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
       br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
       jr(3,4)=n_y*(damp_r*br1+am*br3+dpip_r)
       jr(3,4)=jr(3,4)+t_y*damp_r*br2+am*t_y*br4
c-------------------------------------------------------------
c  ------  f44444444444444444444444444444 ---------
c-------------------------------------------------------------
       hr_l=rold_l*(uold_l*uold_l+vold_l*vold_l)/2.0d0+ga*pold_l/gm1
       hr_r=rold_r*(uold_r*uold_r+vold_r*vold_r)/2.0d0+ga*pold_r/gm1
       f(4)=am*(mpl_m*hr_l+mmin_m*hr_r)
c---------------------
       br1=d2mr_l*hr_l+mpl_m*(uold_l*uold_l+vold_l*vold_l)/2.0d0
       br1=br1+d3mr_l*hr_r
       jl(4,1)=damr_l*f(4)/am+am*br1
c---------------------
       br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
       br1=br1+d3mu_l*hr_r
       jl(4,2)=damu_l*f(4)/am+am*br1
c---------------------
       br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
       br1=br1+d3mv_l*hr_r
       jl(4,3)=damv_l*f(4)/am+am*br1
c---------------------
       br1=d2mp_l*hr_l+mpl_m*ga/gm1
       br1=br1+d3mp_l*hr_r
       jl(4,4)=damp_l*f(4)/am+am*br1
c----------------------------------------------------------
c----------------------------------------------------------
       br1=d2mr_r*hr_l+mmin_m*(uold_r*uold_r+vold_r*vold_r)/2.0d0
       br1=br1+d3mr_r*hr_r
       jr(4,1)=damr_r*f(4)/am+am*br1
c---------------------
       br1=d2mu_r*hr_l+mmin_m*uold_r*rold_r
       br1=br1+d3mu_r*hr_r
       jr(4,2)=damu_r*f(4)/am+am*br1
c---------------------
       br1=d2mv_r*hr_l+mmin_m*vold_r*rold_r
       br1=br1+d3mv_r*hr_r
       jr(4,3)=damv_r*f(4)/am+am*br1
c---------------------
       br1=d2mp_r*hr_l+mmin_m*ga/gm1
       br1=br1+d3mp_r*hr_r
       jr(4,4)=damp_r*f(4)/am+am*br1
c-------------------------------------------------------------
c-------------------------------------------------------------
       return
       end