      PROGRAM GSAPRES                                                                                                               
C     PROGRAM GSAPRES (GSFLD,       GSLNSP,       GSPRES,        GSDPDE,        J2                                                  
C    1                                             INPUT         OUTPUT,)       J2                                                  
C    2            TAPE1=GSFLD, TAPE2=GSLNSP, TAPE3=GSPRES, TAPE4=GSDPDE,                                                            
C    3                                       TAPE5=INPUT,  TAPE6=OUTPUT)                                                            
C     ---------------------------------------------------------------           J2                                                  
C                                                                               J2                                                  
C     APR 07/09 - S.KHARIN: FIX A BUG IN CALCULATION OF DP/DE.                  J2                                                  
C                 (MAXIMUM INACCURACY IN THE OLDER CODE IS OF THE ORDER PLID/P0 J2                                                  
C                  WHICH IS ABOUT 0.1% FOR PLID=50, AND SMALLER FOR HIGHER LID) J2                                                  
C     NOV 17/08 - S.KHARIN                                                                                                          
C                                                                               J2                                                  
CGSAPRES - CALCULATES PRESSURE AND GRADIENT AT SIGMA/HYBRID LEVELS.     2  2 C  J1                                                  
C                                                                               J3                                                  
CAUTHOR  - S.KHARIN                                                             J3                                                  
C                                                                               J3                                                  
CPURPOSE - CALCULATES PRESSURE AND GRADIENT AT SIGMA/HYBRID LEVELS.             J3                                                  
C                                                                               J3                                                  
CINPUT FILES...                                                                 J3                                                  
C                                                                               J3                                                  
C      GSFLD  = SETS OF ETA (SIGMA/HYBRID) LEVEL GRID DATA.                     J3                                                  
C      GSLNSP = SERIES OF GRIDS OF LN(SFC PRES IN MB).                          J3                                                  
C                                                                               J3                                                  
COUTPUT FILE...                                                                 J3                                                  
C                                                                               J3                                                  
C      GSPRES = PRESSURE VALUES ON ETA LEVELS (PA).                             J3                                                  
C      GSDPDE = (OPTIONAL) VERTICAL PRESSURE GRADIENT D PRES/D ETA.             J3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               J5                                                  
C      ICOORD = 4H SIG/ 4H ETA FOR SIGMA/ETA VERTICAL COORDINATES.              J5                                                  
C      PTOIT  = PRESSURE (PA) AT THE LID OF MODEL.                              J5                                                  
C                                                                               J5                                                  
CEXAMPLE OF INPUT CARDS...                                                      J5                                                  
C                                                                               J5                                                  
C*GSAPRES. ET15       50.                                                       J5                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK,LDPDE                                                                                                              
                                                                                                                                    
      INTEGER LEV(100),LEVP(100),KBUF(8)                                                                                            
                                                                                                                                    
      REAL ETA  (100),A    (100),B      (100)                                                                                       
      REAL SIG  (100),FSIG (100),DFLNSIG(101),DLNSIG (100)                                                                          
      REAL PR   (100),PRLOG(100)                                                                                                    
      REAL DA   (100),DB   (100)                                                                                                    
                                                                                                                                    
      REAL PSL(18528),F(1852800)                                                                                                    
      COMMON/ICOM  /IBUF(8),IDAT(18528)                                                                                             
                                                                                                                                    
      DATA MAXX/18528/, MAXL/1852800/, MAXLEV/100/, P0/101320.E0/                                                                   
C---------------------------------------------------------------------                                                              
      NFIL=6                                                                                                                        
      CALL JCLPNT(NFIL,1,2,3,4,5,6)                                                                                                 
      NFIL=NFIL-2                                                                                                                   
      LDPDE=.FALSE.                                                                                                                 
      DO 110 N=1,3                                                                                                                  
  110 REWIND N                                                                                                                      
      IF(NFIL.EQ.4)THEN                                                                                                             
        WRITE(6,'(A)')' ALSO CALCULATE DP/DE.'                                                                                      
        LDPDE=.TRUE.                                                                                                                
        REWIND 4                                                                                                                    
      ENDIF                                                                                                                         
                                                                                                                                    
C     * READ THE CONTROL CARDS.                                                                                                     
                                                                                                                                    
      READ(5,5010,END=900) ICOORD,PTOIT                                         J4                                                  
      IF(ICOORD.EQ.NC4TO8("    ")) ICOORD=NC4TO8(" SIG")                                                                            
      WRITE(6,6010) ICOORD,PTOIT                                                                                                    
                                                                                                                                    
C     * GET ETA VALUES FROM THE GSFLD FILE.                                                                                         
                                                                                                                                    
      CALL FILEV (LEV,NSL,IBUF,1)                                                                                                   
      IF(NSL.LT.1 .OR. NSL.GT.MAXLEV) CALL         XIT('GSAPRES',-1)                                                                
      NWDS = IBUF(5)*IBUF(6)                                                                                                        
      IF(NSL*NWDS.GT.MAXL) CALL                    XIT('GSAPRES',-2)                                                                
      DO 116 I=1,8                                                                                                                  
  116 KBUF(I)=IBUF(I)                                                                                                               
      CALL LVDCODE(ETA,LEV,NSL)                                                                                                     
      DO 118 L=1,NSL                                                                                                                
  118 ETA(L)=ETA(L)*0.001E0                                                                                                         
                                                                                                                                    
C     * EVALUATE THE PARAMETERS OF THE ETA VERTICAL DISCRETIZATION.                                                                 
                                                                                                                                    
      CALL COORDAB (A,B, NSL,ETA, ICOORD,PTOIT)                                                                                     
C                                                                                                                                   
C     * EVALUATE THE PARAMETERS FOR VERTICAL DERIVATIVE D PRES / D ETA                                                              
C                                                                                                                                   
      IF(LDPDE)THEN                                                                                                                 
        IF(ICOORD.EQ.NC4TO8(" SIG") ) THEN                                                                                          
          DO K=1,NSL                                                                                                                
            DA(K)=0.E0                                                                                                              
            DB(K)=1.E0                                                                                                              
          ENDDO                                                                                                                     
        ELSE                                                                                                                        
          IF(ICOORD.EQ.NC4TO8("ET10")) THEN                                                                                         
            POW=1.E0                                                                                                                
          ELSEIF(ICOORD.EQ.NC4TO8(" ETA")) THEN                                                                                     
            POW=2.E0                                                                                                                
          ELSEIF(ICOORD.EQ.NC4TO8("ET15")) THEN                                                                                     
            POW=1.5E0                                                                                                               
          ELSEIF(ICOORD.EQ.NC4TO8("ET16")) THEN                                                                                     
            POW=1.6E0                                                                                                               
          ELSE                                                                                                                      
            CALL                                   XIT('GSAPRES',-3)                                                                
          ENDIF                                                                                                                     
          ETOIT=PTOIT/P0                                                                                                            
          DO K=1,NSL                                                                                                                
            DB(K)=POW/(1.E0-ETOIT)*B(K)**((POW-1.E0)/POW)                                                                           
            DA(K)=P0*(1.E0-DB(K))                                                                                                   
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
                                                                                                                                    
                                                                                                                                    
C---------------------------------------------------------------------                                                              
C     * GET NEXT SET FROM FILE GSFLD.                                                                                               
                                                                                                                                    
      NSETS=0                                                                                                                       
  150 CALL GETSET2 (1,F,LEV,NSL,IBUF,MAXX,OK)                                                                                       
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6030) NSETS                                                                                                         
        IF(NSETS.EQ.0)THEN                                                                                                          
          CALL                                     XIT('GSAPRES',-4)                                                                
        ELSE                                                                                                                        
          CALL                                     XIT('GSAPRES',0)                                                                 
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSAPRES',-5)                                                                
      NAME=IBUF(3)                                                                                                                  
      NPACK=IBUF(8)                                                                                                                 
      NST=IBUF(2)                                                                                                                   
                                                                                                                                    
C     * GET LN(SF PRES) FOR THIS STEP, PUT AT BEGINNING OF COMMON BLOCK.                                                            
                                                                                                                                    
      CALL GETFLD2 (2, PSL ,NC4TO8("GRID"),NST,NC4TO8("LNSP"),1,                                                                    
     +                                           IBUF,MAXX,OK)                                                                      
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
      IF(.NOT.OK) CALL                             XIT('GSAPRES',-6)                                                                
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSAPRES',-7)                                                                
                                                                                                                                    
C     * CALCULATE PRESSURE ON ETA LEVELS.                                                                                           
                                                                                                                                    
      DO K=1,NSL                                                                                                                    
        DO I=1,NWDS                                                                                                                 
          F((K-1)*NWDS+I)=A(K)+B(K)*100.E0*EXP(PSL(I))                                                                              
        ENDDO                                                                                                                       
      ENDDO                                                                                                                         
                                                                                                                                    
C     * WRITE THE PRESSURE LEVEL GRIDS ONTO FILE 3.                                                                                 
                                                                                                                                    
      IBUF(3)=NAME                                                                                                                  
      IBUF(8)=NPACK                                                                                                                 
      CALL PUTSET2 (3,F,LEV,NSL,IBUF,MAXX)                                                                                          
                                                                                                                                    
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
      NSETS=NSETS+1                                                                                                                 
      IF(.NOT.LDPDE)GO TO 150                                                                                                       
                                                                                                                                    
C     * CALCULATE PRESSURE ON ETA LEVELS.                                                                                           
                                                                                                                                    
      DO K=1,NSL                                                                                                                    
        DO I=1,NWDS                                                                                                                 
          F((K-1)*NWDS+I)=DA(K)+DB(K)*100.E0*EXP(PSL(I))                                                                            
        ENDDO                                                                                                                       
      ENDDO                                                                                                                         
                                                                                                                                    
C     * WRITE THE PRESSURE LEVEL GRIDS ONTO FILE 3.                                                                                 
                                                                                                                                    
      IBUF(3)=NAME                                                                                                                  
      IBUF(8)=NPACK                                                                                                                 
      CALL PUTSET2 (4,F,LEV,NSL,IBUF,MAXX)                                                                                          
                                                                                                                                    
      GO TO 150                                                                                                                     
                                                                                                                                    
C     * E.O.F. ON INPUT.                                                                                                            
                                                                                                                                    
  900 CALL                                         XIT('GSAPRES',-8)                                                                
C---------------------------------------------------------------------                                                              
 5010 FORMAT(10X,1X,A4,E10.0)                                                   J4                                                  
 6010 FORMAT(' ICOORD=',1X,A4,', P.LID (PA)=',E10.3)                                                                                
 6030 FORMAT('0 GSAPRES CALCULATED',I5,' SETS OF ',A4)                                                                              
 6035 FORMAT(' ',A4,I10,2X,A4,I10,4I6)                                                                                              
      END                                                                                                                           
