      PROGRAM GSATHLT                                                                                                               
C     PROGRAM GSATHLT (GSTEMP,       GSPHI,       GSRGASM,       GSLNSP,        J2                                                  
C    1                               GTTEMP,       GTPHI,        GTPRES,        J2                                                  
C    2                               GTDPDTH,      INPUT,        OUTPUT,)       J2                                                  
C    3          TAPE11=GSTEMP,TAPE12=GSPHI,TAPE13=GSRGASM,TAPE14=GSLNSP,                                                            
C    4                        TAPE15=GTTEMP,TAPE16=GTPHI, TAPE17=GTPRES,                                                            
C    5                        TAPE18=GTDPDTH,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 (CORRECT TEMP/PHI/PRES/DPDTH REFERENCES                                                                 
C                 IN PUTSET2 CALLS; INTRODUCE COLUMN ARRAYS TEMPSC, PHISC                                                           
C                 INTO CALL TO EATHLT).                                         J2                                                  
C     JUL 27/94 - J. KOSHYK                                                                                                         
C                                                                               J2                                                  
CGSATHLT - INTERPOLATES TEMP, PHI, PRES AND DPRES/DTHETA FROM ETA               J1                                                  
C          (SIGMA/HYBRID) OR PRESSURE LEVELS TO THETA LEVELS.           4  4 C  J1                                                  
C                                                                               J3                                                  
CAUTHOR  - J. KOSHYK                                                            J3                                                  
C                                                                               J3                                                  
CPURPOSE - INTERPOLATES TEMP, PHI, PRES AND DPRES/DTHETA FROM ETA               J3                                                  
C          (SIGMA/HYBRID) OR PRESSURE LEVELS TO THETA LEVELS.                   J3                                                  
C          EXTRAPOLATION UP AND DOWN IS BY LAPSE RATES, -DT/DZ SPECIFIED        J3                                                  
C          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      GSTEMP  = SERIES OF GRIDS OF TEMPERATURE ON ETA OR PRESSURE LEVELS.      J3                                                  
C      GSPHI   = SERIES OF GRIDS OF PHI ON ETA OR PRESSURE LEVELS.              J3                                                  
C      GSRGASM = SERIES OF GRIDS OF MOIST GAS CONSTANT ON ETA OR PRESSURE       J3                                                  
C                LEVELS.                                                        J3                                                  
C      GSLNSP = SERIES OF GRIDS OF LN(SF PRES); NOT USED IF "ICOORD.EQ.4HPRES". J3                                                  
C                                                                               J3                                                  
COUTPUT FILES...                                                                J3                                                  
C                                                                               J3                                                  
C      GTTEMP  = THETA LEVEL TEMPERATURES.                                      J3                                                  
C      GTPHI   = THETA LEVEL GEOPOTENTIALS.                                     J3                                                  
C      GTPRES  = THETA LEVEL PRESSURES.                                         J3                                                  
C      GTDPDTH = THETA LEVEL DPRES/DTHETA.                                      J3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                                                                                   
C      NTHL   = NUMBER OF REQUESTED THETA LEVELS (MAX $L$).                     J5                                                  
C      RLUP   = LAPSE RATE, -DT/DZ USED TO EXTRAPOLATE UPWARDS (DEG/M).         J5                                                  
C      RLDN   = LAPSE RATE USED TO EXTRAPOLATE DOWNWARDS.                       J5                                                  
C      ICOORD = 4H SIG/4H ETA/4HET10/4HET15 FOR INPUT ETA 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*GSATHLT.    5        0.   6.5E-03  ETA        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),PHI(1871328),RGASM(1871328),LNSP(18528)                                                                    
      REAL PRES(1871328),DPDTH(1871328)                                                                                             
                                                                                                                                    
C     * WORKSPACE ARRAYS                                                                                                            
                                                                                                                                    
      REAL TEMPSC(100),PHISC(100)                                                                                                   
      REAL TPRESS(100),TTHS(100),TRMEAN(101),TGAMMA(101),                                                                           
     1     TNU(101)                                                                                                                 
                                                                                                                                    
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
                                                                                                                                    
      DATA MAXX/18528/, MAXL/1871328/, MAXLEV/100/                                                                                  
C---------------------------------------------------------------------------------                                                  
      NFIL=10                                                                                                                       
      CALL JCLPNT(NFIL,11,12,13,14,15,16,17,18,5,6)                                                                                 
      DO 110 N=11,18                                                                                                                
  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('GSATHLT',-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('GSATHLT',-2)                                                                
                                                                                                                                    
C     * GET ETA VALUES FROM THE GSTEMP FILE.                                                                                        
                                                                                                                                    
      CALL FILEV (LEV,NSL,IBUF,11)                                                                                                  
      IF(NSL.LT.1 .OR. NSL.GT.MAXLEV) CALL         XIT('GSATHLT',-3)                                                                
      NWDS = IBUF(5)*IBUF(6)                                                                                                        
      IF((MAX0(NSL,NTHL)+1)*NWDS.GT.MAXL) CALL     XIT('GSATHLT',-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 MULTI-LEVEL TEMPERATURE FIELD FROM FILE GSTEMP.                                                                         
                                                                                                                                    
      NSETS=0                                                                                                                       
  150 CALL GETSET2(11,TEMP,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('GSATHLT',-5)                                                                
        ELSE                                                                                                                        
          CALL                                     XIT('GSATHLT', 0)                                                                
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSATHLT',-6)                                                                
      NPACK = IBUF(8)                                                                                                               
                                                                                                                                    
C     * GET MULTI-LEVEL GEOPOTENTIAL FIELD FROM FILE GSPHI.                                                                         
                                                                                                                                    
      CALL GETSET2(12,PHI,LEV,NSL,IBUF,MAXX,OK)                                                                                     
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
      IF(.NOT. OK) CALL                            XIT('GSATHLT',-7)                                                                
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSATHLT',-8)                                                                
                                                                                                                                    
C     * GET MULTI-LEVEL MOIST GAS CONSTANT FIELD FROM FILE GSRGASM.                                                                 
                                                                                                                                    
      CALL GETSET2(13,RGASM,LEV,NSL,IBUF,MAXX,OK)                                                                                   
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
      IF(.NOT. OK) CALL                            XIT('GSATHLT',-9)                                                                
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSATHLT',-10)                                                               
                                                                                                                                    
C     * GET LN(SF PRES).                                                                                                            
                                                                                                                                    
      NST= IBUF(2)                                                                                                                  
      IF (ICOORD.NE.NC4TO8("PRES")) THEN                                                                                            
       CALL GETFLD2 (14,LNSP,NC4TO8("GRID"),NST,NC4TO8("LNSP"),1,                                                                   
     +                                              IBUF,MAXX,OK)                                                                   
       IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                            
       IF(.NOT. OK) CALL                           XIT('GSATHLT',-11)                                                               
       CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                               
       IF(.NOT.OK) CALL                            XIT('GSATHLT',-12)                                                               
      ELSE                                                                                                                          
       DO 200 I=1,NWDS                                                                                                              
         LNSP(I)=0.0E0                                                                                                              
  200  CONTINUE                                                                                                                     
      ENDIF                                                                                                                         
                                                                                                                                    
C     * INTERPOLATE IN-PLACE FROM ETA TO THETA.                                                                                     
                                                                                                                                    
      CALL EATHLT (TEMP,PHI,RGASM,LNSP,NWDS,NSL,NSL1,TH,                                                                            
     1             NTHL,RLUP,RLDN,A,B,TEMP,PHI,PRES,DPDTH,                                                                          
     2             TEMPSC,PHISC,TPRESS,TTHS,TRMEAN,TGAMMA,TNU)                                                                      
                                                                                                                                    
C     * WRITE THE THETA LEVEL GRIDS ONTO FILES 15, 16, 17 AND 18.                                                                   
                                                                                                                                    
      IBUF(8)=NPACK                                                                                                                 
                                                                                                                                    
      IBUF(3)=NC4TO8("TEMP")                                                                                                        
      CALL PUTSET2 (15, TEMP,LEVTH,NTHL,IBUF,MAXX)                                                                                  
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
                                                                                                                                    
      IBUF(3)=NC4TO8(" PHI")                                                                                                        
      CALL PUTSET2 (16,  PHI,LEVTH,NTHL,IBUF,MAXX)                                                                                  
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
                                                                                                                                    
      IBUF(3)=NC4TO8("PRES")                                                                                                        
      CALL PUTSET2 (17, PRES,LEVTH,NTHL,IBUF,MAXX)                                                                                  
      IF(NSETS.EQ.0) WRITE(6,6035) IBUF                                                                                             
                                                                                                                                    
      IBUF(3)=NC4TO8("DPTH")                                                                                                        
      CALL PUTSET2 (18,DPDTH,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('GSATHLT',-13)                                                               
  909 CALL                                         XIT('GSATHLT',-14)                                                               
C---------------------------------------------------------------------------------                                                  
 5010 FORMAT(10X,I5,2E10.0,1X,A4,E10.0)                                         J4                                                  
 5020 FORMAT(16I5)                                                              J4                                                  
 6010 FORMAT(' RLUP,RLDN = ',2F6.4,', 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 GSATHLT INTERPOLATED',I5,' SETS OF ',A4)                                                                            
 6035 FORMAT(' ',A4,I10,2X,A4,I10,4I6)                                                                                              
      END                                                                                                                           
