Powder toy atmos update code

    for (y=0; y<YRES/CELL; y++) //update velocity and pressure
    {
        for (x=0; x<XRES/CELL; x++)
        {
            dh = 0.0f;
            dx = 0.0f;
            dy = 0.0f;
            for (j=-1; j<2; j++)
            {
                for (i=-1; i<2; i++)
                {
                    if (y+j>0 && y+j<YRES/CELL-2 &&
                            x+i>0 && x+i<XRES/CELL-2 &&
                            !bmap_blockairh[y+j][x+i])
                    {
                        f = kernel[i+1+(j+1)*3];
                        dh += hv[y+j][x+i]*f;
                        dx += vx[y+j][x+i]*f;
                        dy += vy[y+j][x+i]*f;
                    }
                    else
                    {
                        f = kernel[i+1+(j+1)*3];
                        dh += hv[y][x]*f;
                        dx += vx[y][x]*f;
                        dy += vy[y][x]*f;
                    }
                }
            }
            tx = x - dx*0.7f;
            ty = y - dy*0.7f;
            i = (int)tx;
            j = (int)ty;
            tx -= i;
            ty -= j;
            if (i>=2 && i<XRES/CELL-3 && j>=2 && j<YRES/CELL-3)
            {
                odh = dh;
                dh *= 1.0f - AIR_VADV;
                dh += AIR_VADV*(1.0f-tx)*(1.0f-ty)*(bmap_blockairh[j][i] ? odh : hv[j][i]);
                dh += AIR_VADV*tx*(1.0f-ty)*(bmap_blockairh[j][i+1] ? odh : hv[j][i+1]);
                dh += AIR_VADV*(1.0f-tx)*ty*(bmap_blockairh[j+1][i] ? odh : hv[j+1][i]);
                dh += AIR_VADV*tx*ty*(bmap_blockairh[j+1][i+1] ? odh : hv[j+1][i+1]);
            }
            pv[y][x] += (dh-hv[y][x])/5000.0f;
            if(!gravityMode){ //Vertical gravity only for the time being
                float airdiff = hv[y-1][x]-hv[y][x];
                if(airdiff>0 && !bmap_blockairh[y-1][x])
                    vy[y][x] -= airdiff/5000.0f;
            }
            ohv[y][x] = dh;
        }
    }