      PROGRAM LAYERS 
C     PARAMETER (NLAYRS=37,NINTFC=NLAYRS+1,NSTND=12)
      PARAMETER (NLAYRS=38,NINTFC=NLAYRS+1,NSTND=14,LTROP=8)
      PARAMETER (DETA20=.0024285896)
      PARAMETER (T0=288.,G=9.8,R=287.04,GAMMA=.0065
     1,          GORG=G/(R*GAMMA),RGOG=1./GORG)
C     PARAMETER (A=0.8,B=1.03636364,C=-0.97272727,D=0.136363636)
C     PARAMETER (A=0.8,B=1.03636364,C=-1.17272727,D=0.216363636)
C---> PARAMETER (A=0.8,B=1.03636364,C=-1.17272727,D=0.165318)
CCCC      PARAMETER (A=0.8,B=1.03636364,C=-1.17272727,D=0.1647617)
      PARAMETER (A=0.8,B=1.03636364,C=-1.17272727,D=0.1508623)
      REAL EINTFC(NINTFC),DETA(NLAYRS),STNDPR(NSTND),ZETA(NLAYRS)
      INTEGER LDT1(NLAYRS)
      DATA STNDPR/20.,30.,50.,70.,100.,150.,200.,250.,300.,400.
     1,           500.,700.,850.,1000./
      DATA PTOP/25./
C-----------------------------------------------------------------
C***
C***  THIS PROGRAM WILL CALCULATE ETA (SIGMA) LAYER INTERFACES
C***  ACCORDING TO NORM'S ANALYTIC RELATIONSHIP
C***
C-----------------------------------------------------------------
C     READ(5,10)A,B,C,D
   10 FORMAT(4F10.6)
      EITOP=0.
      SUM=0.
      DO 100 N=1,NINTFC
      NI=N-1
      FRAC=REAL(NI)/NLAYRS
      EINTFC(N)=FRAC*(A+FRAC*(B+FRAC*(C+FRAC*D)))
      IF(N.GT.1)THEN
        DETA(N-1)=EINTFC(N)-EITOP
        EITOP=EINTFC(N)
        SUM=SUM+DETA(N-1)
      ENDIF
  100 CONTINUE
      DO 125 N=1,NLAYRS
      DETA(N)=DETA(N)/SUM
  125 CONTINUE
C***
C***  NOW TAKE ENOUGH FROM A TROPOPAUSE LAYER (L=LTROP) TO RAISE 
C***  THE THICKNESS OF THE LOWEST LAYER UP TO 20 METERS
C***
      PIECE=DETA20-DETA(NLAYRS)
      DETA(LTROP)=DETA(LTROP)-PIECE
      DETA(NLAYRS)=DETA(NLAYRS)+PIECE
      DO 130 N=2,NINTFC
      EINTFC(N)=EINTFC(N-1)+DETA(N-1)
  130 CONTINUE
C***
C***  IF THE USE OF 'STANDARD' INPUT PRESSURE LEVELS IS NEEDED,
C***  THE LDT1 ARRAY WILL HOLD THE ETA LAYER INDEX THAT IS
C***  NEAREST TO EACH STANDARD PRESSURE LEVEL.
C***
      DETAMH=0.
      PRES=PTOP
      PD=1013.25-PTOP
      NSTND2=NSTND-2
      NN=1
      DO 200 N=1,NLAYRS
      DETAH=0.5*DETA(N)
      PRES=PRES+(DETAMH+DETAH)*PD
      DETAMH=DETAH
      DIF1=ABS(PRES-STNDPR(NN))
      NN1=MIN(NN+1,NSTND)
      DIF2=ABS(PRES-STNDPR(NN1))
      IF(DIF1.LT.DIF2)THEN
        LDT1(N)=MAX(NN-1,1)
      ELSE
        LDT1(N)=MIN(NN,NSTND2)
        IF(NN.LT.NSTND)NN=NN+1
      ENDIF
  200 CONTINUE
C***
C***  WRITE THE LAYER INFORMATION BOTH FOR VIEWING
C***  AND ALSO FOR INPUT INTO THE NEXT ROUTINE THAT
C***  WILL INSERT ADDITIONAL LAYERS WHEREVER THEY
C***  ARE DESIRED
C***
      SUMETA=0.
      DO 500 N=1,NINTFC
      NI=N-1
      IF(N.EQ.1)THEN
        WRITE(6,450)NI,EINTFC(N)
  450   FORMAT(' INTFC=',I2,'  ETA=',E13.6)
      ELSE
        AETA=EINTFC(NI)+0.5*DETA(NI)
        SUMETA=SUMETA+DETA(NI)
        WRITE(6,460)NI,EINTFC(N),DETA(NI),AETA,LDT1(NI)
  460   FORMAT(' INTFC=',I2,'  ETA=',E13.6,'  DETA=',E13.6,
     1         '  AETA=',E13.6,'  LDT1=',I2)
        WRITE(12,470)DETA(NI)
  470   FORMAT(E13.6)
        WRITE(14,475)NI,DETA(NI),SUMETA
  475   FORMAT(1X,I3,5X,E13.6,5X,E13.6)
      ENDIF
  500 CONTINUE
C***
C***  FIND THE HEIGHTS OF THE ETA INTERFACES FOR A
C***  STANDARD ATMOSPHERE
C***
      DO 550 N=NLAYRS,1,-1
      ZETA(N)=T0*(1.-((PTOP+EINTFC(N)*PD)/(PD+PTOP))**RGOG)/GAMMA
  550 CONTINUE
C
      DO 610 N=1,NLAYRS
      PRS=PD*DETA(N)
      WRITE(6,605)N,PRS,ZETA(N)
  605 FORMAT(' N=',I2,' DELPRS=',E12.5,' ZETA=',E12.5)
  610 CONTINUE
      STOP
      END
