      PROGRAM GSATHL                                                                                                                
C     PROGRAM GSATHL (GSFLD,       GSLNSP,       GSTEMP,       GTHFLD,          J2                                                  
C    1                                            INPUT,       OUTPUT,  )       J2                                                  
C    2          TAPE1=GSFLD, TAPE2=GSLNSP, TAPE3=GSTEMP, TAPE4=GTHFLD,                                                              
C    3                                      TAPE5=INPUT, TAPE6=OUTPUT)                                                              
C     ----------------------------------------------------------------          J2                                                  
C                                                                               J2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       J2                                                  
C     APR 01/98 - F.MAJAESS (REPLACE HOLLERITH IN FORMAT STATEMENTS)                                                                
C     NOV 14/95 - C.READER (DISABLE READING "LNSP" IF "ICOORD.NE.4HPRES")                                                           
C     AUG 16/95 - F.MAJAESS (REPLACE AMAX1 CALLS BY MAX)                                                                            
C     OCT 05/94 - J. KOSHYK (INTRODUCE COLUMN ARRAY FSIG INTO CALL TO EAPL).                                                        
C     JUL 27/94 - J. KOSHYK                                                                                                         
C                                                                               J2                                                  
CGSATHL   - INTERPOLATES A FIELD FROM SIGMA/HYBRID OR PRESSURE LEVELS           J1                                                  
C           TO NTHL THETA (POTENTIAL TEMPERATURE) LEVELS.               3  1 C  J1                                                  
C                                                                               J3                                                  
CAUTHOR  - J. KOSHYK                                                            J3                                                  
C                                                                               J3                                                  
CPURPOSE - INTERPOLATES FROM ETA (SIGMA/HYBRID) OR PRESSURE LEVELS TO NTHL      J3                                                  
C          POTENTIAL TEMPERATURE (THETA) LEVELS. THE INTERPOLATION IS LINEAR    J3                                                  
C          IN LN(THETA).  EXTRAPOLATION UP AND DOWN IS BY LAPSE RATES,          J3                                                  
C          DF/D(LN THETA) SPECIFIED BY THE USER.                                J3                                                  
C          NOTE: AN ARGUMENT FOR GSLNSP FILE HAS TO BE SPECIFIED ON THE         J3                                                  
C                PROGRAM CALL EVEN IF "ICOORD.EQ.4HPRES".                       J3                                                  
C                                                                               J3                                                  
CINPUT FILES...                                                                 J3                                                  
C                                                                               J3                                                  
C      GSFLD  = SETS OF ETA (SIGMA/HYBRID) OR PRESSURE LEVEL GRID DATA.         J3                                                  
C      GSLNSP = SERIES OF GRIDS OF LN(SF PRES); NOT USED IF "ICOORD.EQ.4HPRES". J3                                                  
C      GSTEMP = SERIES OF GRIDS OF TEMPERATURE.                                 J3                                                  
C                                                                               J3                                                  
COUTPUT FILE...                                                                 J3                                                  
C                                                                               J3                                                  
C      GTHFLD = SETS OF THETA LEVEL GRID DATA.                                  J3                                                  
C                                                                                                                                   
C      NTHL   = NUMBER OF REQUESTED THETA LEVELS (MAX $L$).                     J5                                                  
C      RLUP   = LAPSE RATE, (D FLD/D LN(THETA)) USED TO EXTRAPOLATE UPWARDS,    J5                                                  
C      RLDN   = LAPSE RATE USED TO EXTRAPOLATE DOWNWARDS,                       J5                                                  
C      ICOORD = 4H SIG/4H ETA/4HET10/4HET15  FOR INPUT VERTICAL COORDINATES.    J5                                                  
C               4HPRES FOR INPUT PRESSURE COORDINATES.                          J5                                                  
C      PTOIT  = PRESSURE (PA) AT THE LID OF MODEL.                              J5                                                  
C      LEVTH  = THETA LEVELS (K) (MONOTONE DECREASING I.E. TOP OF               J5                                                  
C               ATMOSPHERE TO BOTTOM).                                          J5                                                  
C                                                                               J5                                                  
CEXAMPLE OF INPUT CARDS...                                                      J5                                                  
C                                                                               J5                                                  
C*GSATHL.     5        0.        0.  SIG        0.                              J5                                                  
C*850  700  550  400  330                                                       J5                                                  
C---------------------------------------------------------------------------------                                                  
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
                                                                                                                                    
      INTEGER LEV(100),LEVTH(100),KBUF(8)                                                                                           
                                                                                                                                    
      REAL ETA(100),A(100),B(100),TH(100)                                                                                           
      REAL TEMP(1871328)                                                                                                            
                                                                                                                                    
C     * WORKSPACE ARRAYS                                                                                                            
                                                                                                                                    
      REAL FSIG(100),THS(100),DFDLNTH(101),DLNTH(100)                                                                               
                                                                                                                                    
      COMMON/BLANCK/F(1871328)                                                                                                      
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
                                                                                                                                    
      DATA MAXX/18528/, MAXL/1871328/, MAXLEV/100/                                                                                  
C---------------------------------------------------------------------------------                                                  
      NFIL=6                                                                                                                        
      CALL JCLPNT(NFIL,1,2,3,4,5,6)                                                                                                 
      DO 110 N=1,4                                                                                                                  
  110 REWIND N                                                                                                                      
                                                                                                                                    
C     * READ THE CONTROL CARDS.                                                                                                     
                                                                                                                                    
      READ(5,5010,END=908) NTHL,RLUP,RLDN,ICOORD,PTOIT                          J4                                                  
      IF(ICOORD.EQ.NC4TO8("    ")) ICOORD=NC4TO8(" SIG")                                                                            
      IF(ICOORD.EQ.NC4TO8(" SIG")) THEN                                                                                             
        PTOIT=MAX(PTOIT,0.00E0)                                                                                                     
      ELSE                                                                                                                          
        PTOIT=MAX(PTOIT, 0.100E-09)                                                                                                 
      ENDIF                                                                                                                         
      IF(NTHL.GT.MAXLEV) CALL                      XIT('GSATHL',-1)                                                                 
      READ(5,5020,END=909) (LEVTH(I),I=1,NTHL)                                  J4                                                  
      WRITE(6,5020) (LEVTH(I),I=1,NTHL)                                                                                             
                                                                                                                                    
C     * DECODE LEVELS.                                                                                                              
                                                                                                                                    
      CALL LVDCODE(TH,LEVTH,NTHL)                                                                                                   
                                                                                                                                    
      WRITE(6,6010) RLUP,RLDN,ICOORD,PTOIT                                                                                          
      CALL WRITLEV(TH,NTHL,' TH ')                                                                                                  
                                                                                                                                    
      DO 114 L=2,NTHL                                                                                                               
  114 IF(TH(L).GE.TH(L-1)) CALL                    XIT('GSATHL',-2)                                                                 
                                                                                                                                    
C     * GET ETA VALUES FROM THE GSFLD FILE.                                                                                         
                                                                                                                                    
      CALL FILEV (LEV,NSL,IBUF,1)                                                                                                   
      IF(NSL.LT.1 .OR. NSL.GT.MAXLEV) CALL         XIT('GSATHL',-3)                                                                 
      NWDS = IBUF(5)*IBUF(6)                                                                                                        
      IF((MAX0(NSL,NTHL)+1)*NWDS.GT.MAXL) CALL     XIT('GSATHL',-4)                                                                 
      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.                                                                 
                                                                                                                                    
      IF (ICOORD.NE.NC4TO8("PRES")) THEN                                                                                            
        CALL COORDAB (A,B, NSL,ETA, ICOORD,PTOIT)                                                                                   
      ELSE                                                                                                                          
        DO 119 L=1,NSL                                                                                                              
          A(L) = ETA(L)*100000.E0                                                                                                   
          B(L) = 0.0E0                                                                                                              
  119 CONTINUE                                                                                                                      
      END IF                                                                                                                        
C---------------------------------------------------------------------------------                                                  
C     * GET NEXT SET FROM FILE GSFLD.                                                                                               
                                                                                                                                    
      NSETS=0                                                                                                                       
      N=NWDS+1                                                                                                                      
  150 CALL GETSET2 (1,F(N),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('GSATHL',-5)                                                                 
        ELSE                                                                                                                        
          CALL                                     XIT('GSATHL', 0)                                                                 
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSATHL',-6)                                                                 
      NAME=IBUF(3)                                                                                                                  
      NPACK=IBUF(8)                                                                                                                 
                                                                                                                                    
C     * GET LN(SF PRES) FOR THIS STEP, PUT AT BEGINNING OF COMMON BLOCK.                                                            
                                                                                                                                    
      NST= IBUF(2)                                                                                                                  
      IF (ICOORD.NE.NC4TO8("PRES")) THEN                                                                                            
       CALL GETFLD2 (2, F ,NC4TO8("GRID"),NST,NC4TO8("LNSP"),1,                                                                     
     +                                            IBUF,MAXX,OK)                                                                     
       IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                            
       IF(.NOT.OK) CALL                            XIT('GSATHL',-7)                                                                 
       CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                               
       IF(.NOT.OK) CALL                            XIT('GSATHL',-8)                                                                 
      ELSE                                                                                                                          
       DO 200 I=1,NWDS                                                                                                              
         F(I)=0.0E0                                                                                                                 
  200  CONTINUE                                                                                                                     
      ENDIF                                                                                                                         
                                                                                                                                    
                                                                                                                                    
C     * GET MULTI-LEVEL TEMPERATURE FIELD FROM FILE GSTEMP.                                                                         
                                                                                                                                    
      CALL GETSET2(3,TEMP,LEV,NSL,IBUF,MAXX,OK)                                                                                     
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
      IF(.NOT. OK) CALL                            XIT('GSATHL',-9)                                                                 
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSATHL',-10)                                                                
                                                                                                                                    
C     * INTERPOLATE IN-PLACE FROM ETA TO THETA.                                                                                     
                                                                                                                                    
      CALL EATHL  (F(N),NWDS,TH,NTHL,F(N),NSL,F,RLUP,RLDN,TEMP,                                                                     
     1             A,B,NSL+1,FSIG,THS,DFDLNTH,DLNTH)                                                                                
                                                                                                                                    
C     * WRITE THE THETA LEVEL GRIDS ONTO FILE 4.                                                                                    
                                                                                                                                    
      IBUF(3)=NAME                                                                                                                  
      IBUF(8)=NPACK                                                                                                                 
      CALL PUTSET2 (4,F(N),LEVTH,NTHL,IBUF,MAXX)                                                                                    
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
                                                                                                                                    
      NSETS=NSETS+1                                                                                                                 
      GO TO 150                                                                                                                     
                                                                                                                                    
C     * E.O.F. ON INPUT.                                                                                                            
                                                                                                                                    
  908 CALL                                         XIT('GSATHL',-11)                                                                
  909 CALL                                         XIT('GSATHL',-12)                                                                
C---------------------------------------------------------------------------------                                                  
 5010 FORMAT(10X,I5,2E10.0,1X,A4,E10.0)                                         J4                                                  
 5020 FORMAT(16I5)                                                              J4                                                  
 6010 FORMAT(' RLUP,RLDN = ',2F6.2,', ICOORD=',1X,A4,                                                                               
     1       ', P.LID (PA)=',E10.3)                                                                                                 
 6015 FORMAT('0   ETA LEVELS',I5,15F6.3/(18X,15F6.3))                                                                               
 6020 FORMAT('0 THETA LEVELS',I5,15F6.0/(18X,15F6.0))                                                                               
 6030 FORMAT('0 GSATHL INTERPOLATED',I5,' SETS OF ',A4)                                                                             
 6035 FORMAT(' ',A4,I10,2X,A4,I10,4I6)                                                                                              
      END                                                                                                                           
