      SUBROUTINE DIMSOL (IERROR)

C

C***********************************************************************

C

C   ROUTINE TO SOLVE DUDENEY MODEL EQUATIONS FOR PARAMETERS BETA,HT,FT,A,HF

C   USED IN PLASMA FREQUENCY FUNCTION FN(H), FROM IONOSPHERIC DATA INPUT.

C

C   REFERENCE: J. ATMOS. TERR. PHYS., VOL.40, PAGE 195-203, DUDENEY 1978.

C

C   NUMERICAL SOLUTIONS DEVELOPED AND ROUTINES WRITTEN 10/4/78 BY

C   RICHARD I. KRESSMAN, BRITISH ANTARCTIC SURVEY, MADINGLEY ROAD,

C   CAMBRIDGE  CB3 0ET, ENGLAND.

C

C        INPUT DATA:        (COMMON BLOCK "IONDAT", ALL S.P. VARIABLES)

C               DEFINITION                         NOTES

C   FOF2   F2-LAYER ORDINARY CRITICAL (MHZ)    MANDATORY

C   FOF1   F1-  "       "        "    ( " )    NO F1 GENERATED IF FOF1=<0.0

C   FOE    E -  "       "        "    ( " )    SET TO 0.6 IF =<0.0

C   HFF2   MINIMUM F2 VIRTUAL HEIGHT  (KM)     MANDATORY IF YMF2=<0.0

C   M3000  F2 M(3000) FACTOR                   MANDATORY IF HMF2=<0.0

C   HME    MINIMUM E VIRTUAL HEIGHT   (KM)     SET TO 110. IF =<0.0

C   YME    E-REGION SEMI-THICKNESS    (KM)     SET TO 15. IF =<0.0

C   HMF2   F-REGION PEAK HEIGHT       (KM)     CALCULATED IF =<0.0

C   YMF2   F-REGION SEMI-THICKNESS    (KM)          "      "   "

C

C        OUTPUT DATA:      (COMMON BLOCK "PARMS")

C                  TYPE                            NOTES

C   BETA    )

C   HT      )

C   FT      )      REAL            PLASMA FREQUENCY MODEL PARAMETERS

C   A       )

C   HF      )

C   LCOND          LOGICAL         F1-REGION L-CONDITION EXISTS

C

C   IERROR         INTEGER         ERROR CODE:

C                                   0 - NO ERROR DETECTED

C                                   1 - NO REAL SOLUTION FOR FOF1

C                                   2 - TANGENCY CONDITION VIOLATED

C                                   3 - FOF2<FOE - NO REAL SOLUTION

C                                   4 - VALLEY OUT OF BOUNDS

C

C   REV. 1-12-78  E-F VALLEY ADDED.

C        4- 1-79  D-E EXPONENTIAL ADDED.

C......JAN 81 BY LEO MCNAMARA MADE SINGLE PRECISION FOR CYBER

c......Jul 86 yme changed from 20 to 15

c-----	last modified 3-mar-87

C

C***********************************************************************

C

      LOGICAL LCOND

      REAL M3000

      COMMON/IONDAT/FOF2, FOF1, FOE, HFF2, M3000, HME, YME,

     *              HMF2, YMF2, H0, XD, DV

      COMMON/PARMS/BETA, HT, FT, A, HF, AD, TAU, HD, LCOND

      DATA PI/3.141592654E0/

C

CHECK INPUT DATA, SUBSTITUTE DEFAULTS

C

10    IERROR = 0

      IF (FOE.LE.0.0) FOE = 0.6

      IF (FOF1 .GT. 0.0) DV = 0.0

      IF (DV .LT. 0.0 .OR. DV .GT. 1.0) GO TO 130

      IF (FOF2.LE.FOE) GO TO 120

      IF (HME.LE.0.0) HME = 110.0

      IF (YME.LE.0.0) YME = 15.0

      IF (XD .LE. 0.0) H0 = HME-YME

      IF (HMF2.LE.0.0) HMF2 = HMAX(M3000,FOF2,FOE)

      IF (YMF2.LE.0.0) YMF2 = YMAX(HMF2,FOF2,FOE,HFF2)

C

CHECK MODEL CAN PRODUCE REAL ANSWER

C

20    ZM = (HMF2-HME)/YMF2

      IF (ZM.LT.ACOS(FOE/FOF2)) GO TO 110

C

C  D-E REGION EXPONENTIAL PARAMETERS

C

      HD = H0

      AD = 1E-10

      TAU = 0.0

      IF (XD .EQ. 0.0) GO TO 30

      X = 1.0 - XD*XD

      TAU = XD*XD*XD * YME*FOE/((HME-H0) * SQRT(X) - YME*X)

      HD = HME - YME*SQRT(X)

      AD = (HD-H0) * EXP(TAU/(XD*FOE))

C

CALCULATE BETA,FT,HT USING STAIRCASE ROOT-FINDING ROUTINE

C

C  (REF. ACTON, "NUMERICAL METHODS THAT WORK", PP42-50.)

C

C  SET UP STARTING POINT

30    FT = (FOE+FOF2)/2.0

35    FA = FT

C

C  FIND NEW FT,BETA

C

      BETA = (ZM-ACOS(FT/FOF2))/ACOS(FOE/FT)

      XE = FOE*FOE*(BETA*BETA-1.0)

      XF = BETA*FOE*FOF2

      FT = SQRT((SQRT(XE*XE + 4.0*XF*XF)-XE)/2.0)

c	type*,' fa,ft ',fa,ft

C

CHECK AGAINST LAST ESTIMATE OF FT FOR CONVERGENCE

C

      IF (ABS(FT-FA)/FT.GT.1E-5) GO TO 35

40    HT = HMF2-YMF2*ACOS(FT/FOF2)

C

C  SET UP FOR F1-LAYER PARAMETERS AND TEST FOR F1 PRESENCE

C

50    A = 0.0

      HF = 0.0

      LCOND = .FALSE.

      XE = (HT-HME)/(BETA*YMF2)

      IF (FOF1 .GT. 0.0) GO TO 60

      IF (DV .LE. 0.0) GO TO 90

C

C  VALLEY - EVALUATE A, HF USING FIRST DIFFERENTIAL

C

      XF = 1.0 - DV

      X = 0.3

      QUAD = 0.0

      GO TO 70

C

CALCULATE A,HF USING STAIRCASE ROUTINE

C

60    XF = FOF1/FOE

      XE2 = XE/2.0

      XE3 = 1.0/COS(XE2)

      XP = (XE/(2.0*PI))**2

C

CHECK FOR REAL SOLUTION FOR FOF1

C

      IF (XF.LE.(XE3+2.0*XP*(1.0+SIN(XE2)**2)*XE3**3)) GO TO 100

      X = ACOS(1.0/XF)/XE

      IF (X.GT.1.0) X = 1.0

      X = (X+0.5)/2.0

65    XA = X

      X = -(2.0*TAN(XE*X)**2+1.0)*XP*TAN(PI*X)/(XF*COS(XE*X)-1.0)

      X = (ATAN(1.0/X)/PI+1.0)/2.0

      IF (ABS(X-XA)/X.GT.1E-5) GO TO 65

C

C  TEST FOR L-CONDITION  (D(F)/D(H) >0. AT REDUCED HEIGHT X)

C

      DX = 2.0*PI*(XF-1.0/COS(XE*X))/(XE*TAN(PI*X))

     *     + TAN(XE*X)/COS(XE*X)

      LCOND = DX.GT.0.0

      IF (LCOND) GO TO 80

C

C  RE-EVALUATE TO GET HEIGHT OF MAXIMUM

C

      QUAD = 1.0

70    XA = X

      X = QUAD+ATAN(-2.0*PI*(XF*COS(XE*X)-1.0)/(XE*TAN(XE*X)))/PI

      IF (ABS(X-XA)/X.GT.1E-5) GO TO 70

80    A = (XF-1.0/COS(XE*X))/(1.0-COS(2.0*PI*X))

      HF = HME+X*(HT-HME)

90    RETURN

100   IERROR = 1

      RETURN

110   IERROR = 2

      RETURN

120   IERROR = 3

      RETURN

130   IERROR = 4

      RETURN

      END

      FUNCTION HMAX (M3000,FOF2,FOE)

C***********************************************************************

C

C     EMPIRICAL FUNCTION FOR HMF2

C

C      REFERENCE: B.A.S. SCIENTIFIC REPORT NO.88, DUDENEY 1974.

C

C***********************************************************************

C

      REAL M3000

      XE = FOF2/FOE

      DM = 0.253/(XE-1.215)-0.012

      FM = M3000*SQRT((0.0196*M3000*M3000+1.0) /

     *      (1.2967*M3000*M3000-1.0))

      HMAX = 1490.0*FM/(M3000+DM)-176.0

      RETURN

      END

      FUNCTION YMAX (HMF2,FOF2,FOE,HFF2)

C***********************************************************************

C

C     EMPIRICAL FUNCTION FOR YMF2

C

C      REFERENCE: B.A.S. SCIENTIFIC REPORT NO.88, DUDENEY 1974.

C

C***********************************************************************

C

      XE = FOF2/FOE

      YMAX = HMF2-HFF2+(0.93/(XE-1.23)+0.05)*(HMF2-164.0)

c	set CCIR supp 252 lower limit on ymf2

	if(ymax.lt.35.0) ymax = 35.

      RETURN

      END

