	subroutine mCA1992(r,re)
c	2D meridian density model, Carpenter and Anderson 1992, including
c     (1) plasmasphere
c     (2) trough
c     (3) plasmapause
c     (4) ionosphere
c	MLT range: [0,15]
c	
      implicit none
      integer*4 kc
      parameter (kc=10)

      real*8 re
      complex*16 r(3)

      real*8      dndx(kc),dndy(kc),dndz(kc) 
      common/mod2/dndx,dndy,dndz

      integer*4 ir1,ir2,iw1,iw2,iprint 
      common/io/ir1,ir2,iw1,iw2,iprint  

      integer*4  iq(kc),icmax
      real*8     del(kc),vd(kc),vtp(kc),vtz(kc),bet(kc),den(kc)
      real*8     wp2(kc),mas(kc),omega(kc),teden,tiden
      common/db1/del,vd,vtp,vtz,bet,den,wp2,mas,omega,
     1           teden,tiden,iq,icmax

	real*8 earthr
      parameter (earthr=6370.0d3)
      real*8 pi,ev,emas,hmas,e0,mu0,clight,kb,rad
      parameter (pi=3.14159265358979324d0)
      parameter (ev=1.60210d-19,emas=9.10908d-31,hmas=1.67252d-27)
      parameter (e0=8.85442d-12,mu0=4.0d0*pi*1.0d-7)
      parameter (clight=2.99792501d8)
      parameter (kb=1.38054d-23)
      parameter (rad=pi/180.0d0)

	integer*4 icall
	save icall
	data icall/0/ 

      logical*1 firstcall
      save firstcall
      data firstcall/.true./

c     the following saved for the density model      
c     tMLT : MLT sector
c     Kpmax : for seting up inner edge of Plasmapause (P.P.), Lppi
c     DLpp : the width of P.P.
c     Lpp : the center of P.P.
c     r0: the lower geodistance boundary
c     n0: the density at r=r0
c     ptr: the power of density distribution in the trough
c     ipp: =0, set the center of plasmapause as Lpp
c          =-1, set the inner edge as Lpp
c          =+1, set the outer edge as Lpp
      real*8, save :: tMLT,Kpmax,DLpp,Lpp
      real*8, save :: apsp,bpsp,lnbpsp,atr, ptr
      real*8, save :: rpsp0, npsp0, rtr0,ntr0
      integer*4, save :: ipp
      integer*4, save :: ionmod

      real*8 lcore
      real*8 ripp

      real*8 f1,df1dx,df1dy,df1dz
      real*8 f2,df2dx,df2dy,df2dz
      real*8 l0,dl,t,t2,g,dgdl
      real*8 g1,dg1dx,dg1dy,dg1dz,dg1dl
      real*8 g2,dg2dx,dg2dy,dg2dz,dg2dl
      
      real*8 nli,dnlidx,dnlidy,dnlidz,dnlidr, rliterm
      real*8, save:: rlip, drlip1, drlip2,nlip
      real*8 etao,etahe,etah
c	local variable


	real*8 x,y,z,x2,y2,z2,xy2,rr,r2,l,r3
      real*8 dldx,dldy,dldz
      real*8 nnew,dndxnew,dndynew,dndznew
      real*8 neq,dneqdl
      real*8 npsp,ntr,dnpspdl,dntrdl
      real*8 npsp1,dnpsp1dl,npsp2,dnpsp2dl,dl12

	re=earthr

      if (firstcall) then
c     J. K. Hargreaves The solar-Terrestrial environment 1992
c     P229 Figure 6.12
c     1200 km (r0=1.1881), Ne= 4e10 m-3             
              firstcall=.false.
c     norminal case
c     10, 4.0, 0.5, 1.0, 4d9, 1.0, 5d8, 2
		! ionosphere
	      nlip=1e12
      	rlip=90.0d3/re+1.0d0
      	rlip=350d3/re+1.0d0
      	drlip1=140d3/re
      	drlip2=0.5d6/re
                  
              !read(ir1,*) tMLT,Lpp,DLpp,rpsp0,npsp0,rtr0,ntr0,ripp
		  tMLT=10
		  Lpp=450e-2
		  DLpp=70e-2
		  rpsp0=1.01d0
		  npsp0=4d9
		  rtr0=1.01d0
		  ntr0=5d8
		  ripp=2



              ipp=nint(ripp)

              apsp=10d0**(3.9043+6)
              bpsp=10d0**(-0.3145)
              ptr=-4.5d0
              lnbpsp=dlog(bpsp)

              if(tMLT.lt.6.and.tMLT.ge.0) then ! 0<=MLT<6
                  atr=(5800d0+300d0*tMLT)*1d6
              elseif(tMLT.ge.6.and.tMLT.lt.15) then ! 6<=MLT<15
                  atr=(-800d0+1400d0*tMLT)*1d6
              elseif(tMLT.gt.100.and.tMLT.le.1000) then ! e.g., 450 500, beta=-4.5, -5.0
                  ptr=-tMLT*1d-2
                  atr=(-800d0+1400d0*10d0)*1d6*Lpp**(-4.5) ! Npp@nominal
                  atr=atr*Lpp**(-ptr)
              elseif(tMLT.lt.0) then
                  read(ir1,*) atr,ptr
		  elseif(abs(tMLT-9999).lt.1d-7) then
		  	write(iw1,*) 'enter ionospheric parameters'
		  	read(ir1,*) nlip, rlip, drlip1, drlip2
		  	write(iw1,*) nlip, rlip, drlip1, drlip2
                  atr=(-800d0+1400d0*10d0)*1d6 ! mlt=10
		  else
		      stop
              endif
              write(iw1,101) tMLT,Lpp,DLpp,ipp,npsp0,rpsp0,ntr0,rtr0,ptr
 101          format(1x,'MLT =',f6.1, /
     1               1x, 'Lpp=',f6.1, 1x,'DLpp=',f6.1,1x,'ipp=',i3,/
     2               1x, 'plasmasphere n0=',g14.5,3x,'r0=',g14.5, /     
     2               1x, '      trough n0=',g14.5,3x,'r0=',g14.5, /
     3               1x, 'power of density distr. in trough = ', g14.5)   
      endif

	x=dreal(r(1))
	y=dreal(r(2))
	z=dreal(r(3))
	x2=x*x
      y2=y*y
	z2=z*z
	xy2=x2+y2
	r2=x2+y2+z2
	rr=sqrt(r2)
	l=rr*r2/xy2
      r3=rr*r2
	dldx=3.*rr/xy2*x+r3*(-2.)*x/(xy2*xy2)
	dldy=3.*rr/xy2*y+r3*(-2.)*y/(xy2*xy2)
	dldz=3.*rr/xy2*z



      
c     2D plasmasphere density 
c     equatorial density profile from CA1992 model
c     the latitudinal dependence is controlled by 
c     the constant density npsp0 at a low altitude r=psp0
c     density at lower boundary at constant geocentric distance      
c     set the lower boundary to below the ionosphere, or set r0=1.0d0
c     so that the lower boundary inconsistence is not an issue now. 2010/11/24
c
      npsp=apsp*bpsp**l
      dnpspdl=npsp*lnbpsp

c	Sheeley et al., 2001
c	apsp=1390d6*3d0**4.8d0
c	npsp=apsp/l**4.8d0
c	dnpspdl=npsp/l*(-4.8d0)
	
      call eqmap_mp(rpsp0,npsp0,npsp,dnpspdl,x,y,z,f1,df1dx,df1dy,df1dz)


c     2D trough density 
c     equatorial density profile from CA1992 model
c     the latitudinal dependence is controlled by 
c     the constant density ntr0 at a low altitude r=rtr0
c     
      ntr=atr*l**ptr
      dntrdl=ntr*(ptr/l)
      call eqmap_mp(rtr0,ntr0,ntr,dntrdl,x,y,z,f2,df2dx,df2dy,df2dz)

      l0=Lpp
      dl=DLpp*0.5d0
      t=(l-l0)/dl
      t2=t*t
      g=dexp(-t2)
      dgdl=g*(-2.0d0*t)/dl
      if(t.le.0.0d0)then
            g1=1.0d0
            dg1dl=0.0d0
            g2=g
            dg2dl=dgdl
      else
            g1=g
            dg1dl=dgdl
            g2=1.0d0
            dg2dl=0.0d0
      endif

      dg1dx=dg1dl*dldx
      dg1dy=dg1dl*dldy
      dg1dz=dg1dl*dldz

      dg2dx=dg2dl*dldx
      dg2dy=dg2dl*dldy
      dg2dz=dg2dl*dldz

c ipp flag of choice of plasmapase
c ipp = 0 set the center of pp
c ipp = 1 set the inner edge
c ipp = 2 set the outer edge
c 

      if(ipp.eq.1) then
            g2=1.0d0-g1
            dg2dx=-dg1dx
            dg2dy=-dg1dy
            dg2dz=-dg1dz
      endif

      if(ipp.eq.2) then
            g1=1.0d0-g2
            dg1dx=-dg2dx
            dg1dy=-dg2dy
            dg1dz=-dg2dz
      endif


      nnew=( (f1-f2)*(g1-g2)+f1+f2 )*0.5d0
      dndxnew=( (df1dx-df2dx)*(g1-g2) 
     1			    +(f1-f2)*(dg1dx-dg2dx) +df1dx+df2dx )*0.5d0
      dndynew=( (df1dy-df2dy)*(g1-g2) 
     1			    +(f1-f2)*(dg1dy-dg2dy) +df1dy+df2dy )*0.5d0
      dndznew=( (df1dz-df2dz)*(g1-g2) 
     1			    +(f1-f2)*(dg1dz-dg2dz) +df1dz+df2dz )*0.5d0

c     
c 	write(*,*) nnew, nli
c	write(iw1,*) f1,f2,g1,g2,ipp    




c     ionosphere model density nli, depending on r only
c     with a peak density nlip at r=rlip
c     different falling scale length on the two sides arond the peak
c     nlip: the peak density of lower ionosphere
c     rlip: the peak location
c     drlip1: the scale height of lower edge of density peak
c     drlip2: the scale height of higher edge of density peak


c      nlip=1e12
c      rlip=90.0d3/re+1.0d0
c      rlip=350d3/re+1.0d0
c      drlip1=140d3/re
c      drlip2=0.5d6/re
c	write(*,*) rr,rlip,nlip,drlip1,drlip2
      if(rr.gt.rlip) then
              rliterm=(rr-rlip)/drlip2
              nli=nlip*exp(-rliterm*rliterm)
              dnlidr=nli*(-2.0d0*rliterm/drlip2)
              dnlidx=dnlidr*x/rr
              dnlidy=dnlidr*y/rr
              dnlidz=dnlidr*z/rr
      else
              rliterm=(rr-rlip)/drlip1
              nli=nlip*exp(-rliterm*rliterm)
              dnlidr=nli*(-2.0d0*rliterm/drlip1)
              dnlidr=nli/drlip1
              dnlidx=dnlidr*x/rr
              dnlidy=dnlidr*y/rr
              dnlidz=dnlidr*z/rr
      endif



      

c     set up the first two species
c     e- and H+ in either order    
	      den(1:2)=nnew+nli ! A series
	      ! den(1:2)=nnew+nli+1.0d6 ! B series
            dndx(1:2)=dndxnew+dnlidx
            dndy(1:2)=dndynew+dnlidy
            dndz(1:2)=dndznew+dnlidz
      end

      subroutine eqmap_mp(r0,nlal,neq,dneqdl,x,y,z,
     1 n,dndx,dndy,dndz)
c     to construct density model in the meridian plane
c     use the equatorial density distribution to 
c     infer the distribution off the equator
c     assume fixed density nlal at low altitude r=r0
c     current position: (l, lat)
c     equatorial density: neq, dneqdl
c     current density n, dndx,dndy,dndz
c     Lunjin Chen, 2010/4/15
      implicit none
      real*8 r0,nlal
      real*8 x,y,z
      real*8 x2,y2,z2,xy1,xy2,r2,r3
      real*8 rr,drdx,drdy,drdz
      real*8 l,dldx,dldy,dldz
      real*8 lr,dlrdx,dlrdy,dlrdz

      real*8 neq,dneqdl
      real*8 n,dndx,dndy,dndz
      real*8 alp,dalpdl,dalpdx,dalpdy,dalpdz
      real*8 nle,lr0,lnlr,lnlr0
      real*8 latterm



	x2=x*x
      y2=y*y
	z2=z*z
	xy2=x2+y2
      xy1=sqrt(xy2)
	r2=x2+y2+z2
	rr=sqrt(r2)
      r3=rr*r2

c **** drdxyz ****
	drdx=x/rr
	drdy=y/rr
	drdz=z/rr

c **** L & dLdxyz ****
	l=rr*r2/xy2
	dldx=3.*rr/xy2*x+r3*(-2.)*x/(xy2*xy2)
	dldy=3.*rr/xy2*y+r3*(-2.)*y/(xy2*xy2)
	dldz=3.*rr/xy2*z

c **** Lr & dlrdxyz ****
c L/rr=1/cos^2(lambda)      
      lr=l/rr
      lnlr=dlog(lr)
      dLrdx=dLdx/rr+L*(-1.)/r2*drdx
      dLrdy=dLdy/rr+L*(-1.)/r2*drdy
      dLrdz=dLdz/rr+L*(-1.)/r2*drdz

      nle=nlal/neq
      lr0=l/r0
      lnlr0=dlog(lr0)

      alp=dlog(nle)/dlog(lr0)
      latterm=lr**alp
      n=neq*latterm

      dalpdl=-dneqdl/neq/lnlr0-dlog(nle)/(lnlr0*lnlr0*l)
      dalpdx=dalpdl*dldx
      dalpdy=dalpdl*dldy
      dalpdz=dalpdl*dldz

      dndx=dneqdl*dldx*latterm+
     1      n*(alp/lr*dlrdx+lnlr*dalpdx) 
      dndy=dneqdl*dldy*latterm+
     1      n*(alp/lr*dlrdy+lnlr*dalpdy) 
      dndz=dneqdl*dldz*latterm+
     1      n*(alp/lr*dlrdz+lnlr*dalpdz) 

      end subroutine
