      PROGRAM BETAO                                                                                                                 
C     PROGRAM BETAO (GPHIS,      GPHI,      BETA,      PS,      OUTPUT, )       D2                                                  
C    1         TAPE1=GPHIS,TAPE2=GPHI,TAPE3=BETA,TAPE4=PS,TAPE6=OUTPUT)                                                             
C     -----------------------------------------------------------------         D2                                                  
C                                                                               D2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       D2                                                  
C     JAN 07/93 - E. CHAN   (DECODE LEVELS USING LVDCODE)                                                                           
C     JAN 29/92 - E. CHAN   (CONVERT HOLLERITH LITERALS TO ASCII)                                                                   
C     JUL 04/88 - F.MAJAESS (CHANGE TIME FIELD TO '-1' IN GETFLD2)                                                                  
C     MAR   /87 - G.J.BOER                                                                                                          
C                                                                               D2                                                  
CBETAO   - COMPUTES BETA FOR OBSERVED HEIGHTS                           2  2    D1                                                  
C                                                                               D3                                                  
CAUTHOR  - G.J.BOER                                                             D3                                                  
C                                                                               D3                                                  
CPURPOSE - COMPUTES THE FUNCTION "BETA" CORRESPONDING TO A SET                  D3                                                  
C          OF GEOPOTENTIAL AND SURFACE GEOPOTENTIAL GRIDS.                      D3                                                  
C          FORM OF BETA IS CONSISTENT WITH THE USE OF PROGRAM VPINT.            D3                                                  
C          NOTE - BETA MAY EXCEED 1 IN LOWEST LAYER.                            D3                                                  
C                 MINIMUM NUMBER OF LEVELS IS 2, MAX IS $PL$.                   D3                                                  
C                                                                               D3                                                  
CINPUT FILES...                                                                 D3                                                  
C                                                                               D3                                                  
C      GPHIS  = INPUT GEOPOTENTIAL OF TERRAIN                                   D3                                                  
C      GPHI   = INPUT SETS OF GEOPOTENTIAL ON PRESSURE LEVELS                   D3                                                  
C                                                                               D3                                                  
COUTPUT FILES...                                                                D3                                                  
C                                                                               D3                                                  
C      BETA   = "NEW" BETA AT REQUESTED PRESSURE LEVLES FOR EACH                D3                                                  
C               PRESSURE GRID.                                                  D3                                                  
C      PS     = OUTPUT SERIES OF INFERRED SURFACE PRESSURE                      D3                                                  
C-------------------------------------------------------------------------                                                          
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMMON/BLANCK/PHIS(18528),PHIM(18528),PHI(18528),                                                                             
     1              BETAM(18528),BETA(18528),PS(18528)                                                                              
                                                                                                                                    
      LOGICAL OK,SPEC                                                                                                               
      INTEGER LEV(100)                                                                                                              
      REAL PRH(101),P(101),A(101),DP(101),DA(101)                                                                                   
                                                                                                                                    
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      DATA MAXX/18528/                                                                                                              
      DATA MAXLEV/100/,PHMAX/1013.0E0/                                                                                              
C---------------------------------------------------------------------                                                              
      NFF=5                                                                                                                         
      CALL JCLPNT(NFF,1,2,3,4,6)                                                                                                    
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      REWIND 3                                                                                                                      
      REWIND 4                                                                                                                      
                                                                                                                                    
C     * FIND HOW MANY LEVELS THERE ARE.                                                                                             
                                                                                                                                    
      CALL FILEV(LEV,NLEV,IBUF,2)                                                                                                   
      IF(NLEV.LT.2.OR.NLEV.GT.MAXLEV) CALL         XIT('BETAO',-1)                                                                  
      NAME=IBUF(3)                                                                                                                  
      WRITE(6,6025) IBUF                                                                                                            
      WRITE(6,6005) NAME,NLEV,(LEV(I),I=1,NLEV)                                                                                     
                                                                                                                                    
C     * DETERMINE THE FIELD SIZE AND CHECK FOR SPECTRAL INPUT                                                                       
                                                                                                                                    
      KIND=IBUF(1)                                                                                                                  
      SPEC=(KIND.EQ.NC4TO8("SPEC") .OR. KIND.EQ.NC4TO8("FOUR"))                                                                     
      NWDS=IBUF(5)*IBUF(6)                                                                                                          
      IF(SPEC) CALL                                XIT('BETAO',-2)                                                                  
C                                                                                                                                   
C     * GET PRESSURE, A=LN(PRES), HALF LEVEL PRESSURE,                                                                              
C     * LAYER PRESS THICKNESS DP AND DIFFERENCE DA IN LN(PRESS)                                                                     
C     * *** NB *** CHOICE OF PH(1)=P(1)/2 AND PH(NLEV+1)=1013 MB                                                                    
C     *            IS FOR CONSISTENCY OF NEW BETA AND OLD "DELTA"                                                                   
C     *            WITH VPINT.  LOWEST FULL LEVEL PRESSURE MUST BE                                                                  
C     *            LESS THAN 1013 MB.                                                                                               
C                                                                                                                                   
      NLEVM=NLEV-1                                                                                                                  
      NLEVP=NLEV+1                                                                                                                  
C                                                                                                                                   
C     * DECODE LEVELS.                                                                                                              
C                                                                                                                                   
      CALL LVDCODE(P,LEV,NLEV)                                                                                                      
C                                                                                                                                   
      DO 110 L=1,NLEV                                                                                                               
         A(L)=LOG(P(L))                                                                                                             
  110 CONTINUE                                                                                                                      
C                                                                                                                                   
      IF(P(NLEV).GT.PHMAX) CALL                    XIT('BETAO',-3)                                                                  
      PRH(1)=P(1)*0.5E0                                                                                                             
C                                                                                                                                   
         DO 120 L=2,NLEV                                                                                                            
         PRH(L)=(P(L)+P(L-1))*0.5E0                                                                                                 
         DA(L)=A(L)-A(L-1)                                                                                                          
  120    CONTINUE                                                                                                                   
      PRH(NLEVP)=PHMAX                                                                                                              
      DA(NLEVP)=LOG(PHMAX)-A(NLEV)                                                                                                  
C                                                                                                                                   
      DO 130 L=1,NLEV                                                                                                               
         DP(L)=PRH(L+1)-PRH(L)                                                                                                      
  130 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * GET THE MOUNTAINS.                                                                                                          
C                                                                                                                                   
      CALL GETFLD2(1,PHIS,NC4TO8("GRID"),-1,NC4TO8("PHIS"),1,                                                                       
     +                                          IBUF,MAXX,OK)                                                                       
      IF(.NOT.OK) CALL                             XIT('BETAO',-4)                                                                  
      WRITE(6,6025) IBUF                                                                                                            
      NPACK=MIN(2,IBUF(8))                                                                                                          
                                                                                                                                    
C                                                                                                                                   
C     * CALCULATE NEW BETA                                                                                                          
C                                                                                                                                   
      NSETS=0                                                                                                                       
C                                                                                                                                   
  150 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * SPECIAL PROCESSING FOR FIRST LEVEL. SET FIRST LEVEL BETAM= 1,                                                               
C     * PS=PHMAX AND GET GEOPOTENTIAL PHIM.                                                                                         
C                                                                                                                                   
      DO 200 I=1,NWDS                                                                                                               
         PS(I)=PHMAX                                                                                                                
         BETA(I)=1.0E0                                                                                                              
  200 BETAM(I)=1.0E0                                                                                                                
C                                                                                                                                   
      CALL GETFLD2(2,PHIM,NC4TO8("GRID"),-1,NC4TO8(" PHI"),LEV(1),                                                                  
     +                                               IBUF,MAXX,OK)                                                                  
      IF( .NOT.OK ) THEN                                                                                                            
        IF(NSETS.EQ.0)THEN                                                                                                          
          WRITE(6,6010)                                                                                                             
          CALL                                     XIT('BETAO',-5)                                                                  
        ELSE                                                                                                                        
          WRITE(6,6020) NSETS                                                                                                       
          CALL                                     XIT('BETAO',0)                                                                   
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      IF( NSETS.EQ.0) THEN                                                                                                          
        NPACK=MIN(NPACK,IBUF(8))                                                                                                    
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CHECK THAT FIRST LEVEL IS ABOVE ALL MOUNTAINS.                                                                              
C                                                                                                                                   
         DO 210 I=1,NWDS                                                                                                            
  210    IF(PHIS(I).GT.PHIM(I)) CALL               XIT('BETAO',-6)                                                                  
C                                                                                                                                   
C                                                                                                                                   
C     * PROCESS REMAINING LEVELS                                                                                                    
C                                                                                                                                   
      DO 260 L=2,NLEV                                                                                                               
      CALL GETFLD2(2,PHI,NC4TO8("GRID"),-1,NC4TO8(" PHI"),LEV(L),                                                                   
     +                                              IBUF,MAXX,OK)                                                                   
      IF( .NOT.OK ) THEN                                                                                                            
        WRITE(6,6020) NSETS                                                                                                         
        CALL                                       XIT('BETAO',0)                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
      DO 220 I=1,NWDS                                                                                                               
         BETA(I)=1.0E0                                                                                                              
         IF(BETAM(I).LT.1.0E0) BETA(I)=0.0E0                                                                                        
  220 CONTINUE                                                                                                                      
C                                                                                                                                   
      DO 230 I=1,NWDS                                                                                                               
        IF(PHIS(I).GE.PHI(I).AND.PHIS(I).LT.PHIM(I)) THEN                                                                           
          PS(I)=EXP(A(L-1)+(PHIS(I)-PHIM(I))*DA(L)/(PHI(I)-PHIM(I)))                                                                
            IF(PS(I).LE.PRH(L)) THEN                                                                                                
                BETAM(I)=(PS(I)-PRH(L-1))/DP(L-1)                                                                                   
                BETA(I)=0.0E0                                                                                                       
            ELSE                                                                                                                    
                BETA(I)=(PS(I)-PRH(L))/DP(L)                                                                                        
            ENDIF                                                                                                                   
        ENDIF                                                                                                                       
  230 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * PUT BETAM ON FILE 3.                                                                                                        
C                                                                                                                                   
      IBUF(3)=NC4TO8("BETA")                                                                                                        
      IBUF(8)=NPACK                                                                                                                 
      IBUF(4)=LEV(L-1)                                                                                                              
      CALL PUTFLD2(3,BETAM,IBUF,MAXX)                                                                                               
C                                                                                                                                   
C                                                                                                                                   
C     * SPECIAL PROCESSING FOR LAST HALF LAYER - ASSUME SAME Z-P                                                                    
C     * RELATION AS FOR LAYER JUST ABOVE                                                                                            
C                                                                                                                                   
      IF(L.EQ.NLEV) THEN                                                                                                            
        DO 240 I=1,NWDS                                                                                                             
           IF(PHIS(I).LT.PHI(I)) THEN                                                                                               
               PS(I)=EXP(A(L-1)+(PHIS(I)-PHIM(I))*DA(L)                                                                             
     1                 /(PHI(I)-PHIM(I)))                                                                                           
               BETA(I)=(PS(I)-PRH(NLEV))/DP(NLEV)                                                                                   
           ENDIF                                                                                                                    
  240   CONTINUE                                                                                                                    
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CYCLE THE FIELDS                                                                                                            
C                                                                                                                                   
      DO 250 I=1,NWDS                                                                                                               
         BETAM(I)=BETA(I)                                                                                                           
  250 PHIM(I)=PHI(I)                                                                                                                
C                                                                                                                                   
  260 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * PUT LAST BETA ON FILE 3.                                                                                                    
C                                                                                                                                   
      IBUF(3)=NC4TO8("BETA")                                                                                                        
      IBUF(8)=NPACK                                                                                                                 
      IBUF(4)=LEV(NLEV)                                                                                                             
      CALL PUTFLD2(3,BETAM,IBUF,MAXX)                                                                                               
C                                                                                                                                   
C     * PUT PS ON FILE 4.                                                                                                           
C                                                                                                                                   
      IBUF(3)=NC4TO8("  PS")                                                                                                        
      IBUF(8)=NPACK                                                                                                                 
      IBUF(4)=1                                                                                                                     
      CALL PUTFLD2(4,PS,IBUF,MAXX)                                                                                                  
C                                                                                                                                   
      NSETS=NSETS+1                                                                                                                 
      GO TO 150                                                                                                                     
C---------------------------------------------------------------------                                                              
 6005 FORMAT('0 BETAO ON ',A4,I5,' LEVELS =',20I5/(29X,20I5/))                                                                      
 6010 FORMAT('0.. BETAO INPUT FILE IS EMPTY')                                                                                       
 6020 FORMAT('0',I6,' SETS WERE PROCESSED')                                                                                         
 6025 FORMAT(' ',A4,I10,2X,A4,I10,4I6)                                                                                              
      END                                                                                                                           
