def densm(alt, d0, xm, tz, mn3, zn3, tn3, tgn3, mn2, zn2, tn2, tgn2):
    '''
/*      Calculate Temperature and Density Profiles for lower atmos.  */
'''
    xs = [0.0]*10
    ys = [0.0]*10
    y2out = [0.0]*10
    rgas = 831.4
    #rgas = 831.44621    #maybe make this a global constant?
    densm_tmp=d0
    if (alt>zn2[0]): # pragma: no cover
        if(xm==0.0):
            return tz[0]
        else:
            return d0

    #/* STRATOSPHERE/MESOSPHERE TEMPERATURE */
    if (alt>zn2[mn2-1]):
        z=alt
    else:
        z=zn2[mn2-1];
    mn=mn2;
    z1=zn2[0];
    z2=zn2[mn-1];
    t1=tn2[0];
    t2=tn2[mn-1];
    zg = zeta(z, z1);
    zgdif = zeta(z2, z1);

    #/* set up spline nodes */
    for k in range(mn):
        xs[k]=zeta(zn2[k],z1)/zgdif;
        ys[k]=1.0 / tn2[k];
    yd1=-tgn2[0] / (t1*t1) * zgdif;
    yd2=-tgn2[1] / (t2*t2) * zgdif * (pow(((re[0]+z2)/(re[0]+z1)),2.0));

    #/* calculate spline coefficients */
    spline (xs, ys, mn, yd1, yd2, y2out);   #No need to change this
    x = zg/zgdif;
    y = [0.0]
    splint (xs, ys, y2out, mn, x, y);

    #/* temperature at altitude */
    tz[0] = 1.0 / y[0];
    if (xm!=0.0):
        #/* calaculate stratosphere / mesospehere density */
        glb = gsurf[0] / (pow((1.0 + z1/re[0]),2.0));
        gamm = xm * glb * zgdif / rgas;

        #/* Integrate temperature profile */
        yi = [0.0]
        splini(xs, ys, y2out, mn, x, yi);
        expl=gamm*yi[0];
        if (expl>50.0): # pragma: no cover
            expl=50.0

        #/* Density at altitude */
        densm_tmp = densm_tmp * (t1 / tz[0]) * exp(-expl);

    if (alt>zn3[0]):
        if (xm==0.0):
            return tz[0]
        else:
            return densm_tmp

    #/* troposhere / stratosphere temperature */
    z = alt;
    mn = mn3;
    z1=zn3[0];
    z2=zn3[mn-1];
    t1=tn3[0];
    t2=tn3[mn-1];
    zg=zeta(z,z1);
    zgdif=zeta(z2,z1);



    #/* set up spline nodes */
    for k in range(mn):
        xs[k] = zeta(zn3[k],z1) / zgdif;
        ys[k] = 1.0 / tn3[k];
    
    yd1=-tgn3[0] / (t1*t1) * zgdif;
    yd2=-tgn3[1] / (t2*t2) * zgdif * (pow(((re[0]+z2)/(re[0]+z1)),2.0));

    #/* calculate spline coefficients */
    spline (xs, ys, mn, yd1, yd2, y2out);
    x = zg/zgdif;
    y = [0.0]
    splint (xs, ys, y2out, mn, x, y);

    #/* temperature at altitude */
    tz[0] = 1.0 / y[0];
    if (xm!=0.0):
        #/* calaculate tropospheric / stratosphere density */
        glb = gsurf[0] / (pow((1.0 + z1/re[0]),2.0));
        gamm = xm * glb * zgdif / rgas;

        #/* Integrate temperature profile */
        yi = [0.0]
        splini(xs, ys, y2out, mn, x, yi);
        expl=gamm*yi[0];
        if (expl>50.0): # pragma: no cover
            expl=50.0;

        #/* Density at altitude */
        densm_tmp = densm_tmp * (t1 / tz[0]) * exp(-expl);
    
    if (xm==0.0):
        return tz[0];
    else:
        return densm_tmp;
/*      Calculate Temperature and Density Profiles for lower atmos.  */