function y = toppmv(p,t,q,mass_g,iGasUnitIN)

%% modelled on Scott Hannon's toppmv.f file, in the KLAYERS dist
%% uses the same GUNIT values as would be in header.gunit (of a .rtp file)
%% input   p   = pressure        (mb)
%%         t   = temperature     (kelvin)
%%         q   = gas amount      (whatever units)
%%    mass_g   = molar gas mass  (grams)
%%  iGasUnitIN = units           (from .rtp format)
%% expects input LEVELS info, and changes to LEVELS ppmv

%% output : ppmv (gas unit 10)

%% you might want to hand edit this code, to set which svp routine to use.
%% currently it is set to iSVP = 5, which points to wexsvp, which is straight
%% from the klayers code

%% see also http://cimss.ssec.wisc.edu/itwg/groups/rtwg/Profile_Utility/Units_Conversion.f90.html

iSVP = +1;   %use /asl/matlab/aslutil/svp.m --> svp_hw.m
iSVP = +2;   %use /asl/matlab/aslutil/svp_hw.m
iSVP = +3;   %use /asl/matlab/aslutil/svp_gg.m
iSVP = +4;   %use /asl/matlab/aslutil/svp_flatau.m
iSVP = +5;   %use /home/sergio/MATLABCODE/CONVERT_UNITS/wexsvp.m

iSVP = +4;   %use /asl/matlab/aslutil/svp_flatau.m
iSVP = +5;   %use /home/sergio/MATLABCODE/CONVERT_UNITS/wexsvp.m

mdair = 28.966;  % molecular mass of dry air

if iGasUnitIN == 11        %%%in ppbv
  % PPMV = PPBV*1E-3
  y = q*1e-3;

elseif iGasUnitIN == 12    %%%in vmr (volume mixing ratio)
  % PPMV = VMR*1E+6
  y = q*1e6;

elseif iGasUnitIN == 20    %%% mass mixing ratio in g/kg
  % PPMV = MRgkg*((MDAIR*1E-3)/MASSF)*1E+6
  rjunk =1E+3 * mdair/mass_g;
  y = q * rjunk;

elseif iGasUnitIN == 21    %%% mass mixing ratio in g/g
  % PPMV = MRgg*((MDAIR)/MASSF)*1E+6
  rjunk =1E+6 * mdair/mass_g;
  y = q * rjunk;

elseif iGasUnitIN == 30    %% pp in mb
  % PPMV = (PPmb/PIN)*1E+6
  y = q./p * 1e6;

%C            partial pressure in atm
%             ELSEIF (ICODE .EQ. 31) THEN
%C               PPMV = (PPatm*1013.25/PIN)*1E+6
%                DO IL=1,NLEV
%                   MRIN(IL,IG)=(MRIN(IL,IG)*1013.25/PIN(IL))*1E+6
%                ENDDO

elseif (iGasUnitIN == 40) %% relative humidy in percent ... assumes water!
  % PPMV = (RH%/100)*(SVP/PIN)*1E+6
  svp = dosvp(t,iSVP);
  y = q .* svp./p * 1e6/100;

elseif (iGasUnitIN == 41) %% relative humidity (fraction)
  % PPMV = RH*(SVP/PIN)*1E+6
  svp = dosvp(t,iSVP);
  y = q .* svp./p * 1e6;

else
  fprintf(1,'gas unit = %4i \n',iGasUnitIN);
  error('unknown iGasUnitIN');
end
