      PROGRAM BETA                                                                                                                  
C     PROGRAM BETA (PSFC,       BETFIL,       INPUT,       OUTPUT,      )       D2                                                  
C    1        TAPE1=PSFC, TAPE2=BETFIL, TAPE5=INPUT, TAPE6=OUTPUT)                                                                  
C     ------------------------------------------------------------              D2                                                  
C                                                                               D2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       D2                                                  
C     NOV 06/95 - F.MAJAESS (REVISE DEFAULT PACKING DENSITY VALUE)                                                                  
C     JAN 06/92 - E. CHAN (ADD CALLS TO LVDCODE TO DECODE VERTICAL COORDINATE)                                                      
C     JAN 29/92 - E. CHAN (CONVERT HOLLERITH LITERALS TO ASCII)                                                                     
C     MAR 87 - G.J.BOER.                                                                                                            
C                                                                               D2                                                  
CBETA    - COMPUTES THE FUNCTION BETA FOR ALL PRESSURE LEVELS           1  1 C  D1                                                  
C                                                                               D3                                                  
CAUTHOR  - G.J.BOER                                                             D3                                                  
C                                                                               D3                                                  
CPURPOSE - COMPUTES THE FUNCTION "BETA" FOR A FILE OF SURFACE                   D3                                                  
C          PRESSURE GRIDS. PRESSURE VALUES ARE READ FROM A CARD.                D3                                                  
C          NOTE THAT BETA MAY EXCEED 1 IN LOWEST LAYER. FORM OF                 D3                                                  
C          BETA IS CONSISTENT WITH THE USE OF PROGRAM VPINT.                    D3                                                  
C          MAXIMUM NUMBER OF LEVELS IS $PL$.                                    D3                                                  
C                                                                               D3                                                  
CINPUT FILE...                                                                  D3                                                  
C                                                                               D3                                                  
C      PSFC = INPUT SERIES OF SURFACE PRESSURE GRIDS (MB).                      D3                                                  
C                                                                               D3                                                  
COUTPUT FILE...                                                                 D3                                                  
C                                                                               D3                                                  
C      BETFIL = "NEW" COMPUTED BETA FOR EACH REQUESTED PRESSURE LEVEL;          D3                                                  
C               (ONE SET FOR EACH PRESSURE GRID).                               D3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               D5                                                  
C      NLEV = NUMBER OF PRESSURE LEVELS (<= $PL$).                              D5                                                  
C      PR   = VALUE OF THE PRESSURE IN MILLIBARS                                D5                                                  
C                                                                               D5                                                  
CEXAMPLE OF INPUT CARDS...                                                      D5                                                  
C                                                                               D5                                                  
C*BETA       10                                                                 D5                                                  
C*50. 150. 250. 350. 450. 650. 750. 850. 900. 950.                              D5                                                  
C--------------------------------------------------------------------                                                               
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMMON/BLANCK/PS(18528),BET(18528)                                                                                            
C                                                                                                                                   
      LOGICAL OK                                                                                                                    
      INTEGER LEV(100)                                                                                                              
      REAL PR(100),PH(101)                                                                                                          
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      DATA MAXX/18528/                                                                                                              
      DATA MAXLEV/100/,PHMAX/1013.0E0/                                                                                              
C--------------------------------------------------------------------                                                               
      NFF=4                                                                                                                         
      CALL JCLPNT(NFF,1,2,5,6)                                                                                                      
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
C                                                                                                                                   
C     * READ THE PRESSURE LEVELS FROM CARDS AND DECODE VERTICAL COORDINATE.                                                         
C     * TEST THAT MAX LOWEST LEVEL PRESSURE IS  < 1013 MB.                                                                          
C                                                                                                                                   
      READ(5,5010,END=902) NLEV                                                 D4                                                  
      IF((NLEV.LT.2).OR.(NLEV.GT.MAXLEV)) CALL     XIT('BETA',-1)                                                                   
      READ(5,5012,END=903) (LEV(L),L=1,NLEV)                                    D4                                                  
      CALL LVDCODE(PR,LEV,NLEV)                                                                                                     
C                                                                                                                                   
      IF(PR(NLEV).GT.PHMAX) CALL                   XIT('BETA',-2)                                                                   
C                                                                                                                                   
      NLEVM=NLEV-1                                                                                                                  
      NLEVP=NLEV+1                                                                                                                  
C                                                                                                                                   
C     * HALF LEVELS - LOWEST HALF LEVEL PRESSURE IS SET = PHMAX =1013MB.                                                            
C                                                                                                                                   
      PH(1)=PR(1)*0.5E0                                                                                                             
      DO 110 L=2,NLEV                                                                                                               
      PH(L)=(PR(L)+PR(L-1))*0.5E0                                                                                                   
  110 CONTINUE                                                                                                                      
      PH(NLEVP)=PHMAX                                                                                                               
C                                                                                                                                   
C     * GET THE NEXT SURFACE PRESSURE FIELD INTO PS.                                                                                
C                                                                                                                                   
      NR=0                                                                                                                          
  120 CALL GETFLD2(1,PS,NC4TO8("GRID"),-1,-1,-1,IBUF,MAXX,OK)                                                                       
      IF(.NOT.OK)THEN                                                                                                               
        IF(NR.EQ.0) CALL                           XIT('BETA',-3)                                                                   
        WRITE(6,6010) NR                                                                                                            
        CALL                                       XIT('BETA',0)                                                                    
      ENDIF                                                                                                                         
      IF(NR.EQ.0) THEN                                                                                                              
       WRITE(6,6025) IBUF                                                                                                           
       NPACK=MIN(2,IBUF(8))                                                                                                         
      ENDIF                                                                                                                         
      NLG=IBUF(5)                                                                                                                   
      NLAT=IBUF(6)                                                                                                                  
      NWDS=NLG*NLAT                                                                                                                 
C                                                                                                                                   
C     * COMPUTE BETA FOR EACH PRESSURE LEVEL AND SAVE ON FILE 2.                                                                    
C     * NOTE THAT BETA MAY BE >1 AT LOWEST LEVEL.                                                                                   
C                                                                                                                                   
      DO 150 L=1,NLEV                                                                                                               
C                                                                                                                                   
         DO 130 I=1,NWDS                                                                                                            
         BET(I)=1.E0                                                                                                                
         IF(PS(I).LT.PH(L+1).AND.PS(I).GT.PH(L)) BET(I)=(PS(I)-PH(L))                                                               
     1                                               /(PH(L+1)-PH(L))                                                               
         IF(PS(I).LT.PH(L)) BET(I)=0.0E0                                                                                            
C                                                                                                                                   
C        * SPECIAL PROCESSING FOR LAST HALF LAYER                                                                                   
C                                                                                                                                   
         IF(L.EQ.NLEV.AND.PS(I).GT.PH(NLEVP)) BET(I)=(PS(I)-PH(NLEV))                                                               
     1                                          /(PH(NLEVP)-PH(NLEV))                                                               
  130    CONTINUE                                                                                                                   
C                                                                                                                                   
      IPR=INT(PR(L))                                                                                                                
      CALL SETLAB(IBUF,NC4TO8("GRID"),-1,NC4TO8("BETA"),LEV(L),                                                                     
     +                                        NLG,NLAT,0,NPACK)                                                                     
      CALL PUTFLD2(2,BET,IBUF,MAXX)                                                                                                 
      IF(NR.EQ.0) WRITE(6,6025) IBUF                                                                                                
  150 CONTINUE                                                                                                                      
      NR=NR+1                                                                                                                       
      GO TO 120                                                                                                                     
C                                                                                                                                   
C     * E.O.F. ON INPUT.                                                                                                            
C                                                                                                                                   
  902 CALL                                         XIT('BETA',-4)                                                                   
  903 CALL                                         XIT('BETA',-5)                                                                   
C---------------------------------------------------------------------                                                              
 5010 FORMAT(10X,I5)                                                            D4                                                  
 5012 FORMAT(16I5)                                                              D4                                                  
 6010 FORMAT( '0',I6,'  RECORDS READ')                                                                                              
 6025 FORMAT(' ',A4,I10,2X,A4,I10,4I6)                                                                                              
      END                                                                                                                           
