subroutine MarsTropoCalib_func(mf,P0,T0,elev1,elev2,fc,Tc,dnz,mfe,rng_tro,dop_tro)
	! 
	!=========================================================================
	! gfortran -c MarsTropoCalib_func.f90
	!=========================================================================
	!
	! Purpose: This subroutine calculates the zenithal and slant delays as well 
	! as the Doppler error that are resulting from the radio-signal refraction in
	! the Martian troposphere.
	! The calculation assumes an isothermal profile in Mars atmosphere (T(z)=Tsurf)
	! and is based on the Martian refractivity index from Ho et al 2002:
	!		N = Qd*P/T with Qd = 130.6
	!
	! Two mapping functions are available:
	!	1/ Foelsche and Kirchengast 2002 mapping function:
	! 		m(e)= (Re/H+1)*[cos(arcsin(rcos(e)))-rsin(e)] with r=Re/(Re+H)
	!	2/ simplified mapping function: 1/sin(e)
	!
	! INPUTS:
	!	mf:    Key for mapping function
	!	P0:    Surface pressure in mbar
	!	T0:    Surface temperature in K
	!	elev1: Elevation 1 (start) in radian
	!	elev2: Elevation 2 (end) in radian
	!	fc:    Signal frequency in Hz
	!	Tc:    Doppler time count (time difference between end and start) in sec
	!
	! OUTPUTS:
	!	dnz:	Zenithal delay in meter
	!	mfe:	mapping factor 
	!	rng_tro: one-way range error or slant delay in meter
	!	dop_tro: one-way Doppler contribution in Hz
	!
	!	A log file is also generated
	!
	!Historic:
	!---------
	! author: S. Le Maistre
	! Datum: March 2019
	! Last update: SLM - April 1st, 2019 - Correct quantity for rng_tro (was not
	!                                      the slant delay before but the range 
	!                                      difference between slant delay at end
	!                                      and slant delay at start time/elev.)
	!              SLM - October 25, 2019- mean mapping factor is used to compute
	!                                      the slant delay (rng_tro)
	!==========================================================================
	! Variable declaration
	!---------------------

	implicit none

	character(LEN=100)	:: fileOut		! output file
	character*10		:: mf

	double precision	:: P0,fc,Tc,T0
	double precision	:: C,g,M,Rg,rr,Z,Q,H
	double precision	:: elev1,elev2
	double precision	:: dnz,rng_tro
	double precision	:: mfe,mfe1,mfe2,Re
	double precision	:: dop_tro
	double precision	:: Pi, degrad

	intrinsic asin 
	intrinsic acos

	!==========================================================================
	! Constantes
	!-----------

	 Pi=acos(-1.d0)
	 degrad=Pi/180.d0

	 C=299792458.d0   	! ligth speed 
	 Rg=8.31446d0     	! constante des gaz parfait
	 Re=3.396d6       	! Mars equatorial radius
	 g=3.711d0        	! Mars gravitational constant
	 !M=44.01d-3       	! Molar mass of CO2
     M=0.043346d0       ! Molar mass of Mars atmosphere M=0.9532*44.01e-3+0.027*28.02e-3+0.016*39.95e-3
	 Q=130.6d-6       	! Martian refractivity indexe after Ho et al (2002)
	 H=11000.d0         ! Mars troposphere height from Ho et al (~H=Rg*T/(g*M))
     H=H*2.d0           ! = Hatm in Foelsche and Kirchengast 2002
	 z=2.5d5          	! receiver altitude set to above the ionopause

	 rr=Re/(Re+H)
	 
	!======================================================
	! inoutput files
	! ---------------
	 fileOut='MarsTropoCalib_funt.log'
	 open(unit=2,file=fileOut,RECL=2000,position='append')

	!======================================================
  
	 !zenithal delay of Martian atmosphere
	 dnz = Q * Rg * P0 / (g*M) * (1.d0 - exp(-g*M*z/(Rg*T0)))    ! in m

	 !mapping function #1 at t2 - tc/2
	 if (elev1 .ge. 0.d0) then
	    if (mf=="FK2002") then                                           ! FOELSCHE
	        mfe1 = (1.d0+Re/H) * (cos(asin(rr*cos(elev1))) - rr*sin(elev1))
	    else                                                             ! simplified 
	        mfe1 = 1.d0 / sin(elev1)
	    endif
	 else
	    mfe1=0.d0
	 endif

	 !mapping function #2 at t2 + tc/2
	 if (elev2 .ge. 0.d0) then
	    if (mf=="FK2002") then                                           ! FOELSCHE
	        mfe2 = (1.d0 + Re/H) * (cos(asin(rr*cos(elev2))) - rr*sin(elev2))
	    else                                                             ! simplified	 
	        mfe2 = 1.d0 / sin(elev2)
	    endif
	 else
	    mfe2 = 0.d0
	 endif
	 mfe=(mfe1+mfe2)/2.d0	! mean mapping factor

	 !1-way Range error or slant delay due to Mars troposphere along the path 
	 rng_tro = mfe * dnz

	 !1-way Doppler shift due to Mars troposphere along the path
	 dop_tro = fc * (mfe2 - mfe1) * dnz / (C * Tc)
 
	 !writing log file
	 write(2,*)P0,T0,elev1,elev2,dnz,mfe1,mfe2,rng_tro,fc,Tc,dop_tro

end subroutine MarsTropoCalib_func