subroutine vremap_conserve_mass_mixing(src_pressure, src_values, dst_pressure, dst_values)
real, intent(in) :: src_pressure(:,:,:)
real, intent(in) :: src_values(:,:,:)
real, intent(in) :: dst_pressure(:,:,:)
real, intent(inout) :: dst_values(:,:,:)
real, allocatable :: temp_pressures_src(:,:,:), temp_values_src(:,:,:)
real, allocatable :: temp_pressures_dst(:,:,:), temp_values_dst(:,:,:)
real :: src_max_p
integer :: lb_src, lb_dst, lm_src, lm_dst, ub_src, ub_dst, im, jm
lm_src = size(src_values,3)
lm_dst = size(dst_values,3)
lb_src = lbound(src_pressure,3)
lb_dst = lbound(dst_pressure,3)
ub_src = ubound(src_pressure,3)
ub_dst = ubound(dst_pressure,3)
im = size(src_values,1)
jm = size(src_values,2)
! src gets extra level that is zero becasue gmap persists src value in dst below surface
src_max_p = maxval(src_pressure(:,:,ub_src))
allocate(temp_pressures_src(im,jm,lb_src:ub_src+1))
allocate(temp_values_src(im,jm,lm_src+1))
temp_pressures_src(:,:,lb_src:ub_src) = src_pressure
temp_values_src(:,:,1:lm_src) = src_values
temp_pressures_src(:,:,ub_src+1) = src_pressure(:,:,ub_src)+10.0
temp_values_src(:,:,lm_src+1) = 0.0
! add an extra level on dst because if src is below destination we will need the extra stuff
! we need to make sure "extra" stuf from src gets included
allocate(temp_pressures_dst(im,jm,lb_dst:ub_dst+1))
allocate(temp_values_dst(im,jm,lm_dst+1))
temp_pressures_dst(:,:,lb_dst:ub_dst) = dst_pressure
temp_pressures_dst(:,:,ub_dst+1) = src_max_p + 10.0
call gmap(im, jm, lm_src+1, temp_pressures_src, temp_values_src, lm_dst+1, temp_pressures_dst, temp_values_dst)
! add back the "extra" level, have to convert to emission to do it
temp_values_dst(:,:,lm_dst) = temp_values_dst(:,:,lm_dst)*(temp_pressures_dst(:,:,lm_dst+1)-temp_pressures_dst(:,:,lm_dst))/MAPL_GRAV &
+ temp_values_dst(:,:,lm_dst+1)*(temp_pressures_dst(:,:,lm_dst+2)-temp_pressures_dst(:,:,lm_dst+1))/MAPL_GRAV
! add convert back to mass mixing
temp_values_dst(:,:,lm_dst) = temp_values_dst(:,:,lm_dst)*MAPL_GRAV/(temp_pressures_dst(:,:,lm_dst+1)-temp_pressures_dst(:,:,lm_dst))
dst_values = temp_values_dst(:,:,1:lm_dst)
end subroutine vremap_conserve_mass_mixing