C**********************************************************                     
	real function PLAS1(ALT)
C T.L. Gulyaeva....................................................May 2014
CC
CC==== IZMIRAN PLASMASPHERE ELECTRON DENSITY PROFILE ABOVE h05top
CC     Ref. Gulyaeva T.L., X.Huang, and B.W. Reinisch. Plasmaspheric extension of 
CC     topside electron density profiles. Adv. Space Res., 29, No.6, 825-831, 2002. 
CC
CC     T.Gulyaeva............................., July 2004 - September 2013.
CC Formulation is changed fitting plasmasphere model to topside ionosphere model at
CC the Topside Transition Ionosphere Height
CC HEITOP at, Electron Density ELTOP = 0.5*NmF2.
CC Output : electron density, PLAS, [m^-3], at altitude X  [km] above HEITOP
CC 
C
C  Parameters used by function PLAS1:
C
C  W = R12 = Rzs 12-months smoothed sunspot number
C  GLON - geographic longitude
C  HOUR - local time, UT - universal time, 
C  NEW= GMLT Geomagnetic Local Time 
C  LDAYNR - day-of-year 
C  ALON, BLON - corrected geomagnetic coordinates
C  HEITOP = h05top - altitute in the topside at Ne05 = 0.5* NmF2
C  X - altitude in the range:  [h05top, 36,000 km]
C
      INTEGER  ND
      real N6370,L,LT
	real eltop
	COMMON /bplas/w,glat,glon,put,ldaynr,heitop,eltop
      PARAMETER ( PI=3.14159265, RE=6370., RAD=PI/180.)
C
C Specify corrected geomagnetic coordinates alon, blon
C SMI SUBR. ABLON TO PREPARE FOR PLASMASPHERE MODEL
      alon=365.
      blon=365.
      if (((glon.ge.130.).and.(glon.le.178.)).or.((glon.ge.330.).
     &and.(glon.le.340.))) goto 209
      CALL ABLON(GLAT,GLON,ALON,BLON)
  209 continue

	 HHUP=HEITOP
      ND=LDAYNR
	LT=PUT+GLON/15.0
	if (LT.ge.24.) LT=LT-24.0
	if (LT.lt.0.) LT=LT+24.0
      L= 1.0 + ALT/RE
c
c-	A1 = 1.-0.1*COS(PI*(LT-4.)/12.)

C A1-cor - corrected A1(L,GMLT)
	xn = 0.9
	xd = 1.-0.1*cos(PI/L**5)
c-	A1 = xn+(xd-xn)/(1.+exp(7.-GMLT))+(xn-xd)/(1.+exp(17.-GMLT))
	A1 = xn+(xd-xn)/(1.+exp(11.-LT))+(xn-xd)/(1.+exp(21.-LT))  ! LT function

      A2 = 1.+0.2*SQRT(W)
      A3 = 1.+0.001*W
      AL = 1.-EXP(-0.04*L**4)
      AL2 = 1.-EXP(-0.04*2**4)
C !!!!
cold	c1=3.0E+09
C !!New!!    CORRECTION: C1=3E+09 AS UPPER LIMIT OF COEFFICIENT
c-	C1=.0002*ELTOP*(RE-HHUP)

C New expression for C1(GMLT)
	c1=(1.0+cos(-pi*(LT-12.)/24.)/2.)*1.0E+09
C!!!!
      if (alon.le.360.) goto 10
      A0 = (1.0+0.7*COS(RAD*(GLON+70.)))*(1.+COS(2*PI*(ND+16)/365.25))
      BL = 0.5*A0*AL-0.7*L*EXP(0.1*A0*AL)
      BL2 = 0.5*A0*AL2-1.4*EXP(0.1*A0*AL2)
c-      N6370 = 3.0E9*A1*A2*EXP(A3*BL2)
      N6370 = c1*A1*A2*EXP(A3*BL2)
c+!!! old formula:      H = 1342.5/ALOG(CN1000/N6370) IS REPLACED BY :
	H=(RE-HHUP)*.25/ALOG(ELTOP/N6370)
      if (ALT.le.6370.) then
C!!! old formula:    DN = CN1000*EXP((1000.-ALT)/H/L/L) IS REPLACED BY :
	   DN=ELTOP*EXP((HHUP-ALT)/H/L/L)
      else
C!!! old formula:    DN = 3.0E9*A1*A2*EXP(A3*BL) IS REPLACED BY :
	   DN=C1*A1*A2*EXP(A3*BL)
      endif
      Plas1=DN
      RETURN
   10 a0a= (1.0+0.7*COS(RAD*(ALON+70.)))*(1.+COS(2*PI*(ND+16)/365.25))
      b0b= (1.0+0.7*COS(RAD*(BLON+70.)))*(1.+COS(2*PI*(ND+16)/365.25))
      bla= 0.5*a0a*AL-0.7*L*EXP(0.1*a0a*AL)
      blb= 0.5*b0b*AL-0.7*L*EXP(0.1*b0b*AL)
      bl2a= 0.5*a0a*AL2-0.7*2*EXP(0.1*a0a*AL2)
      bl2b= 0.5*b0b*AL2-0.7*2*EXP(0.1*b0b*AL2)
      drea = 3.0E9*A1*A2*EXP(A3*bl2a)
      dreb = 3.0E9*A1*A2*EXP(A3*bl2b)
      drea=alog10(drea)
      dreb=alog10(dreb)
      delon=blon-alon
      if (delon.le.0.) delon=delon+360.
      deglon=glon-alon
      if (deglon.le.0.) deglon=deglon+360.
      dne=drea+(dreb-drea)*deglon/delon
      N6370=10.**dne
C!!!      H = 1342.5/ALOG(CN1000/N6370)
	H=(RE-HHUP)*.25/ALOG(ELTOP/N6370)
      if (ALT.le.6370.) then
C!!!         DN = CN1000*EXP((1000.-ALT)/H/L/L)
	   DN=ELTOP*EXP((HHUP-ALT)/H/L/L)
      else
C!!!         dna = 3.0E9*A1*A2*EXP(A3*bla)
	         dna = C1*A1*A2*EXP(A3*bla)
         dna=alog10(dna)
C!!!         dnb = 3.0E9*A1*A2*EXP(A3*blb)
         dnb = C1*A1*A2*EXP(A3*blb)
         dnb=alog10(dnb)
         dnx=dna+(dnb-dna)*deglon/delon
         dn=10.**dnx
      endif
      Plas1=DN
      RETURN
      END
C
C------------------------------------------------------------------
      SUBROUTINE ablon(xlati,xlongi,alonn,blonn)
C==== GRID POINTS FOR THE FUNCTION PLAS(alt,...) 
      DIMENSION xlat(37),xlon(37),zlon(6)
      DATA xlat /10.62,10.43,10.06,9.60,9.10,8.76,8.60,8.74,9.29,9.82,
     &9.83,9.27,8.38,7.48,7.07,7.05,6.25,4.07,1.76,0.36,-0.42,-1.37,
     &-2.84,-4.01,-4.15,-4.37,-6.34,-9.73,-12.75,-14.46,-14.12,-10.41,
     &-3.90,2.66,7.44,9.95,10.62/
      DATA xlon /0.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.,110.,120.,
     &130.,140.,150.,160.,170.,180.,190.,200.,210.,220.,230.,240.,250.,
     &260.,270.,280.,290.,300.,310.,320.,330.,340.,350.,360./
      do 1 i=1,36
      k=i
      if ((xlongi.ge.xlon(i)).and.(xlongi.lt.xlon(i+1))) goto 2
    1 CONTINUE
    2 fimeq=xlat(k)+(xlat(k+1)-xlat(k))*(xlongi-xlon(k))*0.1
      if ((xlongi.lt.178.).or.(xlongi.gt.330.)) goto 3
      if (xlati.le.(fimeq+1.)) goto 11
      if (xlati.lt.2.6) goto 4
      alonn=178.
      blonn=330.
      goto 11
    4 nl=0
      do 5 i=1,36
c      k=i
      if (((xlati.ge.xlat(i)).and.(xlati.lt.xlat(i+1))).or.
     &((xlati.le.xlat(i)).and.(xlati.gt.xlat(i+1)))) goto 15
      goto 5
   15 nl=nl+1
      zlon(nl)=xlon(i)+10.*(xlati-xlat(i))/(xlat(i+1)-xlat(i))
    5 continue
      alonn=zlon(1)
      blonn=zlon(2)
      goto 11
    3 if (xlati.ge.(fimeq-1.)) goto 11
      if (xlati.gt.7.6) goto 6
      alonn=340.
      blonn=130.
      goto 11
    6 nl=0
      do 7 i=1,36
c      k=i
      if (((xlati.ge.xlat(i)).and.(xlati.lt.xlat(i+1))).or.
     &((xlati.le.xlat(i)).and.(xlati.gt.xlat(i+1)))) goto 17
      goto 7
  17  nl=nl+1
      zlon(nl)=xlon(i)+10.*(xlati-xlat(i))/(xlat(i+1)-xlat(i))
    7 continue
      if ((xlongi.ge.zlon(1)).and.(xlongi.le.130.)) goto 8
      alonn=zlon(nl)
      blonn=zlon(1)
      goto 11
    8 continue
      do 9 i=1,nl
      k=i
      if ((xlongi.ge.zlon(i)).and.(xlongi.le.zlon(i+1))) goto 10
    9 continue
   10 alonn=zlon(k)
      blonn=zlon(k+1)
   11 RETURN
      END
C
C
C**********************************************************
C New:	 T.L. Gulyaeva.................................Sep 30, 2013
	subroutine prepPlas(hm,aep,bb)
C Determine height h05top(0.5*NmF2) for NeQuick Plasmasphere starting height:
      real NeNeQ
	dimension hm(3),aep(3),bb(6)
	COMMON /bplas/r12i,glat,ylon,put,ldaynr,heitop,eltop
C
C Start iterations to obtain h05top:
c	eltop=aep(1)*1E+11/8.   ! 0.5*NmF2
	htop=hm(1)      ! start from hmF2
	eltop=NeNeQ(htop,hm,aep,bb)/2. ! 0.5NmF2
	dh=100.
	iflag=0
    1  htop=htop+dh
    	e1=NeNeQ(htop,hm,aep,bb)
	if (iflag.gt.3) goto 2
	if (e1.lt.eltop) then
	htop=htop-dh
	dh=dh/10.
	iflag=iflag+1
	goto 1
	                  else
	goto 1
	endif
    2	heitop=htop
      RETURN
	end
C***********************************************************************
      real function NeNeQ(h,hm,aep,bb)
      real NeMdGR
      dimension hm(3),aep(3),bb(6)
	COMMON /bplas/r12i,glat,ylon,put,ldaynr,heitop,eltop
 	 if ((h.gt.heitop).and.(heitop.gt.hm(1))) then
C choose option Plasmasphere:
	NeNeq=plas1(h)
	return
	endif
      if (h.gt.hm(1)) then
C choose option topside Ionophere 
         aNmax=NeMdGR(aep,hm,bb,hm(1))
         NeNeQ=topq(h,aNmax,hm(1),bb(6))
         return
      endif
C choose option bottomside Ionophere 
      NeNeQ=NeMdGR(aep,hm,bb,h)
      return
      end
C------------------------------------------------------
      real function NeMdGR(aep,hm,bb,h)
      dimension aep(3),hm(3),bb(6),B(3)

      save

      parameter (f1=10.0, f2=1.0)
      parameter (h0=90.0, Hd=5.0)
      B(1)=bb(5)
      B(2)=bb(3)
      B(3)=bb(1)
      if (h.gt.hm(3)) B(3)=bb(2)
      if (h.gt.hm(2)) B(2)=bb(4)
         sum=0.0
         do jj=1,3
            arg0=(h-hm(jj))
            arg=arg0/B(jj)
            if (jj.gt.1) then
               d=abs(h-hm(1))
               arg=arg*exp(f1/(1.0+f2*d))
            endif
          if (h.lt.h0) arg=arg*(Hd+h0-h)/Hd
            if (abs(arg).gt.25.0) then
               s0=0.0
            else
               ee=exp(arg)
               s0=aep(jj)*ee/(1.0+ee)**2
            endif
            sum=sum+s0
         enddo
         NeMdGR=sum*1.0E11
         return
      end
C-------------------------------------------------------------
      real function topq(h,aNo,hmax,Ho)
      parameter (g=0.125,rfac=100.0)
         dh=h-hmax
         gh=g*dh
         z=dh/(Ho*(1.0+rfac*gh/(rfac*Ho+gh)))
         ee=fexp(z)
         if (ee.gt.1.0E11) then
            ep=4.0/ee
         else
            ep=4.0*ee/(1.0+ee)**2
         endif
         topq=aNo*ep
      return
      end
C-------------------------------------------------------------
	REAL FUNCTION DINTEC(h1,h2,hm,aep,bb)
C==== TECion+TECpl INTEGRATION PROCEDURE TEC/m^3
	REAL NeNeQ
      dimension hm(3),aep(3),bb(6)
	common /bplas/r12i,glat,ylon,gmlt,ldaynr,heitop,eltop
      dnemin=NeNeQ(h1,hm,aep,bb)
c      k=ifix((ah-hmin)/20.+0.5)
C 
	if ((h2-h1).gt.1000.) then
	 st=20.	!  step=20 km above heitop
	else
      st=10.0
	endif
	k=int((h2-h1)/st+0.5)
      x1=h1
      y1=dnemin
      SUM=0.0
      do 1 l=1,k
      x2=x1+st
      y2=NeNeQ(x2,hm,aep,bb)
      SUM=SUM+(y2+y1)*st
	x1=x2
      y1=y2
    1 continue
	dintec=sum*500.
      return
      end
C
C-------------------------------------------------------------
      subroutine subinteg(h2,tec1,hm,aep,bb)
C  integration along vertical profile
C h2 = upper boundary for vTEC [80.0,hpl]
      dimension hm(3),aep(3),bb(6)
		COMMON /bplas/r12i,glat,ylon,put,ldaynr,heitop,eltop
	h0=80.0
	h1=hm(1)
crem	h2=20200.0
        if (h2.le.heitop) then
          tec1=dintec(h0,h2,hm,aep,bb)
	goto 25
        else
	h1a=hm(1)
	h1b=heitop
                tec1=dintec(h0, h1a,hm,aep,bb)
                tec2=dintec(h1a,h1b,hm,aep,bb)
                tec3=dintec(h1b,h2, hm,aep,bb)
                tec4=tec1+tec2+tec3
              endif
  25	tec1=tec4/1.0E15 !TECU/10. equals to 10^15 m^-2
	RETURN
	END
C----------------------------------------------------------------------------
