      SUBROUTINE F2LAYR(IFLAG,FOF2,HMAX,FOE)
C***********************************************************************
C* SUBROUTINE NAME :  F2LAYR            CPC:  SS/MOD/RTN/F2LAYR        *
C*          VERSION:  88-II            DATE:  16 MAR 88                *
C*                                                                     *
C*          PURPOSE:  THIS SUBROUTINE IS THE ELKINS-RUSH MODIFICATION  *
C*                    TO THE URSI MODEL.  GIVEN A SSN AND A KP/QE AS   *
C*                    INPUT, IT WILL CALCULATE THE POSITION OF THE     *
C*                    POLEWARD TROUGH AND MAKE APPROPRIATE MODIFIC-    *
C*                    ATIONS OF THE TROUGH AND AURORAL FOF2 VALUES.    *
C*                                                                     *
C*          CALLING                                                    *
C*        ARGUMENTS:  IFLAG,FOF2,HMAX,FOE                              *
C*                                                                     *
C*    ABBREVIATIONS:  NONE                                             *
C*                                                                     *
C*    COMMON BLOCKS:  F6PARM, GWC1                                     *
C*                                                                     *
C*   FILES ACCESSED:  NONE                                             *
C*                                                                     *
C*           METHOD:  1.  CALCULATE THE LOCAL TIME AND THE POSITION    *
C*                        OF THE TROUGH WALL.                          *
C*                    2.  DETERMINE THE CORRECTION FACTORS FOR THE     *
C*                        POLAR CAVITY,AURORAL ZONE, AND AURORAL TROUGH*
C*                    3.  CALCULATE FOF2 AND HMAX AT THE GRIDPOINT     *
C*                                                                     *
C*          REMARKS:  NONE                                             *
C*                                                                     *
C*       REFERENCES:  1.  ICED-III SYSTEM DOCUMENTATION                *
C*                    2.  SPERRY FORTRAN-77 (ASCII) PRM                *
C*                    3.  AFGWC FORM 10 DX-50194                       *
C*                    4.  AFGWC FORM 10 DX-60096                       *
C*                                                                     *
C* GLOBAL VARIABLES:                                                   *
C*  COMMON /F6PARM/IYR,IMO,IDA,TIME,DAYJUL,EFFSSN,EFFQ,EFFKP,          *
C*                 GGLAT,GGLON,CGLAT,CGLON,TCGM                        *
C*     IYR    IS  YEAR OF DATA SET                                     *
C*     IMO    IS  MONTH OF DATA SET                                    *
C*     IDA    IS  DAY OF DATA SET                                      *
C*     TCGM   RS  CORRECTED GEOMAGNETIC TIME OF A GRIDPOINT            *
C*     TIME   RS  UNIVERSAL TIME IN HH.HH OF DATA SET                  *
C*     DAYJUL RS  JULIAN DATE OF DATA SET                              *
C*     EFFSSN RS  EFFECTIVE SUNSPOT NUMBER                             *
C*     EFFQ   RS  EFFECTIVE AURORAL Q VALUE                            *
C*     EFFKP  RS  EFFECTIVE PLANETARY K INDEX                          *
C*     GGLAT  RS  GEOGRAPHIC LATITUDE OF A GRIDPOINT                   *
C*     GGLON  RS  GEOGRAPHIC LONGITUDE OF A GRIDPOINT                  *
C*     CGLAT  RS  CORRECTED GEOMAGNETIC LATITUDE OF A GRIDPOINT        *
C*     CGLON  RS  CORRECTED GEOMAGNETIC LONGITUDE OF A GRIDPOINT       *
C*                                                                     *
C*  COMMON /GWC1/ KF(10),UF(13,76),KX(10),UX(13,76)                    *
C*     KF     IV  HEADER RECORD VALUES FOR FOF2 COEFFICIENTS IN URSI   *
C*     UF     RA  FOURIER SERIES COEFFICIENTS FOR FOF2 IN URSI         *
C*     KX     IV  HEADER RECORD VALUES FOR M3000 COEFFICIENTS IN URSI  *
C*     UX     RA  FOURIER SERIES COEFFICIENTS FOR M3000 IN URSI        *
C*                                                                     *
C*  LOCAL VARIABLES:  IFLAG    IS  QUALIFIER AT GRIDPOINT              *
C*                    IERR     IS  ERROR FLAG FROM COORDINATE TRANSFORM*
C*                    CHI      RS  SOLAR ZENITH ANGLE                  *
C*                    DELTAN   RS  SMALL INCREMENT OF TIME             *
C*                    DELTAT   RS  MLT HOURS FROM 0300 MLT             *
C*                    EB       RS  EQUATORIAL BOUNDARY VALUE           *
C*                    DUM      RS  PLACEHOLDER IN CALL TO SOLPOS       *
C*                    FOE      RS  FREQUENCY OF E LAYER                *
C*                    FOF      RS  TEMP VALUE OF FREQUENCY OF F2 LAYER *
C*                    FOF2     RS  FREQUENCY OF THE FOF2 LAYER         *
C*                    HMAX     RS  HEIGHT OF THE F2 LAYER              *
C*                    PB       RS  POLAR BOUNDARY OF AURORAL ZONE      *
C*                    PHIA     RS  AURORAL ZONE EQUATORWARD BOUNDARY   *
C*                    PI       RS  3.1415926 - CONSTANT                *
C*                    RDUM     RS  UNUSED RETURN VALUE FROM COORDINATE *
C*                                 TRANSFORMATION                      *
C*                    T        RS  TEMP TIME VARIABLE                  *
C*                    T1       RS  TEMP TIME VARIABLE                  *
C*                    TIMLOC   RS  LOCAL TIME AT GGLON                 *
C*                    TCGM     RS  CORRECTED GEOMAGNETIC TIME          *
C*                    X1       RS  TEMP VARIABLE                       *
C*                    XA       RS  TEMP VARIABLE                       *
C*                    TRFPB    RS  POLEWARE BOUNDARY OF THE TROUGH     *
C*                    TRFMAX   RS  CENTER OF THE TROUGH                *
C*                    TRFEB    RS  EQUATORWARD BOUNDARY OF THE TROUGH  *
C*                    HFRM20   RS  MLT HOURS PAST 20 MLT               *
C*                    CGLN20   RS  CORRECTED GEOMAG LONG. AT 20 MLT    *
C*                    CGLN03   RS  CORRECTED GEOMAG LONG. AT 03 MLT    *
C*                    CGLN07   RS  CORRECTED GEOMAG LONG. AT 07 MLT    *
C*                    GLAT     RS  TEMPORARY GEOGRAPHIC LATITUDE       *
C*                    GLON     RS  TEMPORARY GEOGRAPHIC LONGITUDE      *
C*                    M3000    RS  TEMPORARY M(3000) FACTOR            *
C*                    HMAXTL   RS  F2 HEIGHT AT TOP LEFT OF            *
C*                                 INTERPOLATION BOX                   *
C*                    HMAXBL   RS  F2 HEIGHT AT BOTTOM LEFT OF         *
C*                                 INTERPOLATION BOX                   *
C*                    HMAXTR   RS  F2 HEIGHT AT TOP RIGHT OF           *
C*                                 INTERPOLATION BOX                   *
C*                    HMAXBR   RS  F2 HEIGHT AT BOTTOM RIGHT OF        *
C*                                 INTERPOLATION BOX                   *
C*                    FRLAT    RS  FRACTIONAL LAT. FROM BOTTOM TO TOP  *
C*                                 OF INTERPOLATION BOX                *
C*                    FRLON    RS  FRACTIONAL LONG. FROM LEFT TO RIGHT *
C*                                 OF INTERPOLATION BOX                *
C*                    PT1      RS  TEMP. BOTTOM F2 HEIGHT IN LONGITUDE *
C*                                 INTERPOLATION BOX                   *
C*                    PT2      RS  TEMP. TOP F2 HEIGHT IN LONGITUDE    *
C*                                 INTERPOLATION BOX                   *
C*                    F        RS  M3000 CALCULATION VARIABLE          *
C*                    XE       RS  RATIO OF FOF2 TO FOE                *
C*                    DELTAM   RS  M3000 CORRECTION TERM               *
C*  (I=INTEGER, R=REAL, C=CHARACTER, S=SCALAR, V=VECTOR, A=ARRAY)      *
C*                                                                     *
C*     SETC OPTIONS:  NONE USED                                        *
C*                                                                     *
C*       SUBPROGRAMS                                                   *
C*           CALLED:  QOVAL1 - CALCULATES THE AURORAL OVAL POSITION    *
C*                    FLDHR  - RETRIEVES FOURIER ANALYSIS COEFFICIENTS *
C*                    EQBOUN - CALCULATES THE EQ. BOUNDARY OF OVAL     *
C*                    CNVCRD - CONVERTS GEOMAG TO GEOGRAPHIC COORDS    *
C*                    SOLPOS - COMPUTES SOLAR ZENITH ANGLE             *
C*                                                                     *
C*       SUBPROGRAM                                                    *
C*       CALLED BY :  FIELD6                                           *
C*                                                                     *
C* PROGRAM WRITTEN:  JAN 86 - NATIONAL GEOPHYSICAL DATA CENTER         *
C*                   MAY 86 - AFGWC/SDDE - UPGRADED DOCUMENTATION      *
C*                   AUG 87 - NATIONAL GEOPHYSICAL DATA CENTER (87-I)  *
C*                   MAR 88 - NGDC (88-II) NEW SOLAR ZENITH ANGLE      *
C*                   JAN 89 - NGDC (89)                                *
C*                   APR 90 - NGDC (90-REF)                            *
C*                                                                     *
C***********************************************************************
      INTEGER IDA
      INTEGER IMO
      INTEGER IYR
      REAL CGLAT
      REAL CGLON
      REAL DAYJUL
      REAL EFFKP
      REAL EFFQ
      REAL EFFSSN
      REAL GGLAT
      REAL GGLON
      REAL TCGM
      REAL TIME
      COMMON/F6PARM/IYR,IMO,IDA,TIME,DAYJUL,EFFSSN,EFFQ,EFFKP,
     +              GGLAT,GGLON,CGLAT,CGLON,TCGM
      INTEGER KF
      INTEGER KX
      REAL UF
      REAL UX
      COMMON/GWC1/KF(10),UF(13,76),KX(10),UX(13,76)
 
      INTEGER IFLAG
      INTEGER IERR
      REAL CGLN20
      REAL CGLN03
      REAL CGLN07
      REAL CHI
      REAL DELTAM
      REAL DELTAN
      REAL DELTAT
      REAL EB
      REAL DUM
      REAL F
      REAL FOE
      REAL FOF
      REAL FOF2
      REAL FRLAT
      REAL FRLON
      REAL GLAT
      REAL GLON
      REAL HFRM20
      REAL HMAX
      REAL HMAXBL
      REAL HMAXBR
      REAL HMAXTL
      REAL HMAXTR
      REAL M3000
      REAL PB
      REAL PHIA
      REAL PI
      REAL PT1
      REAL PT2
      REAL RDUM
      REAL T
      REAL T1
      REAL TEMPDM
      REAL TIMLOC
      REAL TRFEB
      REAL TRFMAX
      REAL TRFPB
      REAL X1
      REAL XA
      REAL XE
      DATA PI/3.141593/
 
C.......................................................................
C     AURORAL ZONE BOUNDARIES (CENTER OF AZ = EQ BOUNDARY OF TROUGH)   .
C.......................................................................
      CALL QOVAL1(EFFQ,TCGM,PB,EB,PHIA)
 
C.......................................................................
C     DETERMINE THE CORRECTION FACTOR IN AURORAL ZONE AND POLAR CAVITY .
C.......................................................................
      DELTAT = ABS(TCGM-3.)
      IF (TCGM.GT.15.0) DELTAT = 27. - TCGM
 
      IF (CGLAT.LT.PHIA) GO TO 40
      XA = PB-PHIA
      X1 = (PHIA-CGLAT)/XA
      IF (CGLAT.GT.EB) X1 = X1 * 2.0
      DELTAN=EXP((-1.0*X1*X1)/2.0)
      GO TO 80
C.......................................................................
C     CALCULATE THE CORRECTION FACTOR SOUTH OF AURORAL ZONE            .
C.......................................................................
40    T1 = 0.0
      XA = 3.7 + (1.3*EFFKP)
C.......................................................................
C     COMPUTE SOLAR ZENITH ANGLE                                       .
C.......................................................................
      TIMLOC = TIME + GGLON/15.0
      IF ( TIMLOC.GE.24.00 ) TIMLOC = TIMLOC - 24.00
      IF ( TIMLOC.LT. 0.00 ) TIMLOC = TIMLOC + 24.00
      CALL SOLPOS( IYR,IMO,IDA,TIMLOC,GGLAT,GGLON,CHI,DUM,DUM,DUM )
 
      X1=(CGLAT-PHIA)/XA
      IF (CHI.LE.90.0.OR.(TCGM.GT.6.0.AND.TCGM.LT.18.0))  T1 = 0.0
      IF ((CHI.GT.94.6).AND.((18.0.LE.TCGM).OR.(TCGM.LE.6.0))) T1 =-0.2
      IF ( ((CHI.GE.90.0).AND.(CHI.LE.94.6)).AND.
     +     ((18.0.GE.TCGM).OR.(TCGM.LE.6.0)) ) T1=-0.2 *(CHI-90.0)/4.6
 
C.....................................................................
C     MILLER AND GIBBS 1975                                          .
C.....................................................................
      T=T1*EXP((X1-(X1*X1))/2.0)*EXP((-1.0*DELTAT*DELTAT)/12.0)
      DELTAN=T*(1.0+COS(PI*(DAYJUL+11.0)/182.5))
80    CALL FLDHR(KF,UF,GGLAT,GGLON,FOF,TIME)
      FOF2=FOF*(1.0+DELTAN)
C......................................................................
C     'DENSITY' CORRECTION IN AURORAL ZONE AND POLAR CAVITY CONSIDERED.
C     CONSTANT IN TIME. CORRECTION BASED ON DENSITY DIFFERENCE BETWEEN.
C     4 AND 5 MHZ. 5 MHZ BEING APPROX 1.22663 TIMES 4 MHZ WHICH IS A  .
C     TYPICAL NIGHTTIME URSI FOF2 VALUE.                              .
C......................................................................
      IF (CGLAT.GE.PHIA) FOF2 = SQRT( (FOF*FOF) + (9.0+DELTAN) )
 
C......................................................................
C     HMF2 EQUATION FROM :                                            .
C     DUDENY, J.R., "THE ACCURACY OF SIMPLE METHODS FOR DETERMINING   .
C       THE HEIGHT OF THE MAXIMUM ELECTRON CONCENTRATION OF THE F2-   .
C       LAYER FROM SCALED IONOSPHERIC CHARACTERISTICS",               .
C     J. ATMOS. TERR. PHYSICS, VOL 45, NO. 8/9, PP. 629-640, 1983     .
C     - ERROR IN HPF2 20% WHEN XE=2.0; 60% WHEN XE=1.7                .
C     - MOST APPROPRIATE FOR 15-60 DEGREES MAGNETIC DIP LATITUDE      .
C     - BY CCIR RECOMMENDATION, ADJUST XE TO BE >= 1.7                .
C......................................................................
 
C......................................................................
C     IFLAG >= 100 IMPLIES GRIDPOINT IN THE SUB-AURORAL TROUGH        .
C......................................................................
      IF ( IFLAG.LT.100 ) GO TO 170
 
C......................................................................
C     F2 TROUGH HEIGHT CALCULATION BY INTERPOLATION.                  .
C     INTERPOLATION CASES :                                           .
C       |-------------------------|-----------------------| TRFPB     .
C       |                         |                       |           .
C       |            1            |           3           | ~PHIA     .
C       |                         |                       |           .
C       |-------------------------*-----------------------| TRFMAX    .
C       |                         |                       |           .
C       |            2            |           4           |           .
C       |                         |                       |           .
C       |-------------------------|-----------------------| TRFEB     .
C     20 MLT                    03 MLT                  07 MLT        .
C                                                                     .
C       * = 450.0 KILOMETERS (MAX HEIGHT IN TROUGH)                   .
C......................................................................
      TRFMAX = PHIA - 1.5
      TRFPB = PHIA + 1.5
      CALL EQBOUN( TCGM,TRFMAX,DAYJUL,EFFKP,TRFEB )
      IF ( TRFEB.GT.TRFMAX-1.5 ) TRFEB = TRFMAX - 1.5
 
C......................................................................
C     HOURS PAST 20 MLT                                               .
C......................................................................
      HFRM20 = 4.0 + TCGM
      IF ( TCGM.GE.20.0 ) HFRM20 = TCGM - 20.0
 
C.......................................................................
C     CG LONGITUDES OF 20,03,07 MLT                                    .
C.......................................................................
      CGLN20 = CGLON - ( HFRM20 * 15.0 )
      CGLN03 = CGLN20 + 105.0
      CGLN07 = CGLN20 + 165.0
      IF ( CGLN20.LE.  0.0 ) CGLN20 = CGLN20 + 360.0
      IF ( CGLN03.GE.360.0 ) CGLN03 = CGLN03 - 360.0
      IF ( CGLN07.GE.360.0 ) CGLN07 = CGLN07 - 360.0
 
      XE = FOF2 / FOE
      IF ( XE.LT.1.7 ) XE = 1.7
      DELTAM = ( 0.253/(XE-1.215) ) - 0.012
 
C......................................................................
C     INTERPOLATION CASE 1                                            .
C......................................................................
      IF ( (CGLAT.GE.TRFMAX).AND.(HFRM20.LE.7.) ) THEN
        TEMPDM = DELTAM
        IF ( CGLAT.GT.PHIA ) XE     = AMAX1( XE    ,2.2 )
        DELTAM = ( 0.253/(XE-1.215) ) - 0.012
        CALL CNVCRD( TRFPB,CGLN20,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXTL = (1490.*F)/(M3000+DELTAM) - 176.
        CALL CNVCRD( TRFMAX,CGLN20,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXBL = (1490.*F)/(M3000+TEMPDM) - 176.
        CALL CNVCRD( TRFPB,CGLN03,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXTR = (1490.*F)/(M3000+DELTAM) - 176.
        HMAXBR = 450.0
        FRLAT = (CGLAT-TRFMAX) / (TRFPB-TRFMAX)
        FRLON = HFRM20 / (27.-20.)
      ENDIF
C......................................................................
C     INTERPOLATION CASE 2                                            .
C......................................................................
      IF ( (CGLAT.LT.TRFMAX).AND.(HFRM20.LE.7.) ) THEN
        CALL CNVCRD( TRFMAX,CGLN20,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXTL = (1490.*F)/(M3000+DELTAM) - 176.
        CALL CNVCRD( TRFEB,CGLN20,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXBL = (1490.*F)/(M3000+DELTAM) - 176.
        HMAXTR = 450.0
        CALL CNVCRD( TRFEB,CGLN03,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXBR = (1490.*F)/(M3000+DELTAM) - 176.
        FRLAT = (CGLAT-TRFEB) / (TRFMAX-TRFEB)
        FRLON = HFRM20 / (27.-20.)
      ENDIF
C......................................................................
C     INTERPOLATION CASE 3                                            .
C......................................................................
      IF ( (CGLAT.GE.TRFMAX).AND.(HFRM20.GT.7.) ) THEN
        TEMPDM = DELTAM
        IF ( CGLAT.GT.PHIA ) XE     = AMAX1( XE    ,2.2 )
        DELTAM = ( 0.253/(XE-1.215) ) - 0.012
        CALL CNVCRD( TRFPB,CGLN03,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXTL = (1490.*F)/(M3000+DELTAM) - 176.
        HMAXBL = 450.0
        CALL CNVCRD( TRFPB,CGLN07,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXTR = (1490.*F)/(M3000+DELTAM) - 176.
        CALL CNVCRD( TRFMAX,CGLN07,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXBR = (1490.*F)/(M3000+TEMPDM) - 176.
        FRLAT = (CGLAT-TRFMAX) / (TRFPB-TRFMAX)
        FRLON = (HFRM20-7.) / (31.-27.)
      ENDIF
C......................................................................
C     INTERPOLATION CASE 4                                            .
C......................................................................
      IF ( (CGLAT.LT.TRFMAX).AND.(HFRM20.GT.7.) ) THEN
        HMAXTL = 450.0
        CALL CNVCRD( TRFEB,CGLN03,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXBL = (1490.*F)/(M3000+DELTAM) - 176.
        CALL CNVCRD( TRFMAX,CGLN07,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXTR = (1490.*F)/(M3000+DELTAM) - 176.
        CALL CNVCRD( TRFEB,CGLN07,0.,4,GLAT,GLON,RDUM,2,IERR)
        CALL FLDHR( KX,UX,GLAT,GLON,M3000,TIME )
        F = M3000*SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
        HMAXBR = (1490.*F)/(M3000+DELTAM) - 176.
        FRLAT = (CGLAT-TRFEB) / (TRFMAX-TRFEB)
        FRLON = (HFRM20-7.) / (31.-27.)
      ENDIF
 
C......................................................................
C     TWO-DIMENSIONAL LINEAR INTERPOLATION BETWEEN FOUR CORNERS       .
C......................................................................
      PT1 = HMAXBL + FRLON*(HMAXBR-HMAXBL)
      PT2 = HMAXTL + FRLON*(HMAXTR-HMAXTL)
      HMAX = PT1 + FRLAT*(PT2-PT1)
      RETURN
 
C......................................................................
C     NON-TROUGH HMF2 COMPUTATION                                     .
C......................................................................
170   XE = FOF2 / FOE
      IF ( XE.LT.1.7 ) XE = 1.7
      IF ( (CGLAT.GE.PHIA).AND.(CGLAT.LE.PB) ) XE = AMAX1( XE,2.2 )
      DELTAM = (0.253/(XE-1.215)) - 0.012
      CALL FLDHR(KX,UX,GGLAT,GGLON,M3000,TIME)
      F = M3000 * SQRT((0.0196*M3000*M3000+1)/(1.2967*M3000*M3000-1))
      HMAX = (1490.*F)/(M3000+DELTAM) - 176.
      RETURN
      END
 
