      PROGRAM GSAPZL                                                                                                                
C     PROGRAM GSAPZL (GSPHI,       GSPHIS,       GSLNSP,       GPPHI,           J2                                                  
C    1                             GPTEMP,       INPUT,        OUTPUT,  )       J2                                                  
C    2         TAPE11=GSPHI,TAPE12=GSPHIS,TAPE13=GSLNSP,TAPE14=GPPHI,                                                               
C    3                      TAPE15=GPTEMP,TAPE5 =INPUT, TAPE6 =OUTPUT)                                                              
C     ----------------------------------------------------------------          J2                                                  
C                                                                               J2                                                  
C     MAR 04/10 - S.KHARIN,F.MAJAESS (CHANGE "RGAS" FROM 287. TO 287.04)        J2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                           
C     JAN 12/93 - E. CHAN   (DECODE LEVELS IN 8-WORD LABEL)                                                                         
C     JAN 29/92 - E. CHAN   (CONVERT HOLLERITH LITERALS TO ASCII)                                                                   
C     MAR 14/89 - F.MAJAESS (FIX 6010 FORMAT)                                                                                       
C     FEB 29/88 - R.LAPRISE.                                                                                                        
C                                                                               J2                                                  
CGSAPZL  - CALCULATES T & PHI ON PRESSURE LEVELS FROM PHI ON ETA                J1                                                  
C          (SIGMA/HYBRID) LEVELS                                        3  2 C  J1                                                  
C                                                                               J3                                                  
CAUTHOR  - R. LAPRISE                                                           J3                                                  
C                                                                               J3                                                  
CPURPOSE - CALCULATES T AND PHI ON PRESSURE LEVELS FROM "PHI" ON                J3                                                  
C          ETA (SIGMA/HYBRID) LEVELS.                                           J3                                                  
C          NOTE - THE FOLLOWINGS ARE THE NECESSARY STEPS TO BE TAKEN            J3                                                  
C                 IF TEMP IS USED INSTEAD OF PHI AS INPUT:                      J3                                                  
C                                                                               J3                                                  
C                 1) INPUT FILE... GETEMP, INSTEAD OF GSPHI.                    J3                                                  
C                 2) READ-IN LAYERING PARAMETER LAY.                            J3                                                  
C                 3) FILEV RETURNS SH, TEMPERATURE MID LAYER POSITION,          J3                                                  
C                 4) BASCAL COMPUTES SHB (LAYER BASES) FROM SH AND LAY.         J3                                                  
C                 5) COORDAB COMPUTES PARAMETERS A AND B (ETA COORD)            J3                                                  
C                    FROM SHB.                                                  J3                                                  
C                 6) MODIFY ETAPHI TO EXTRAPOLATE UPWARD USING UPPERMOST        J3                                                  
C                    TEMPERATURE.                                               J3                                                  
C                                                                               J3                                                  
CINPUT FILES...                                                                 J3                                                  
C                                                                               J3                                                  
C      GSPHI  = SETS OF "PHI" ON ETA LEVELS, DEFINED AS                         J3                                                  
C               "PHI"=PHIS + INT(R*T*D LN ETA)  IN PAKGCMH.                     J3                                                  
C      GSPHIS = MOUNTAINS (SURFACE GEOPOTENTIAL).                               J3                                                  
C      GSLNSP = SET OF LN(SURFACE PRESSURE).                                    J3                                                  
C                                                                               J3                                                  
COUTPUT FILES...                                                                J3                                                  
C                                                                               J3                                                  
C      GPPHI  = PRESSURE LEVEL GEOPOTENTIALS.                                   J3                                                  
C      GPTEMP = PRESSURE LEVEL TEMPERATURES.                                    J3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               J5                                                  
C      NPL    = NUMBER OF PRESSURE LEVELS, (MAX $L$).                           J5                                                  
C      RLUP   = DT/D LN(SIG) FOR TEMP EXTRAP UPWARDS.                           J5                                                  
C      GAMMA  = DT/DZ        FOR TEMP EXTRAP DOWNWARD (DEG/M).                  J5                                                  
C      ICOORD = 4H SIG/4H ETA FOR SIGMA/ETA COORDINATES.                        J5                                                  
C      PTOIT  = PRESSURE (PA) OF LID AT TOP OF MODEL.                           J5                                                  
C      PR     = PRESSURE LEVELS (MB)                                            J5                                                  
C                                                                               J5                                                  
CEXAMPLE OF INPUT CARDS...                                                      J5                                                  
C                                                                               J5                                                  
C*GSAPZL.     5        0.    6.5E-3  SIG        0.                              J5                                                  
C*100  300  500  850 1000                                                       J5                                                  
C-------------------------------------------------------------------------                                                          
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
                                                                                                                                    
      INTEGER LG(101),LP(101)                                                                                                       
                                                                                                                                    
      REAL AG   (101),BG   (101),AH    (100),BH     (100)                                                                           
      REAL SG   (101),SH   (100),PR    (100),PRLOG  (100)                                                                           
      REAL SIG  (101),FSIG (101),DLNSIG(101),DFLNSIG(102)                                                                           
                                                                                                                                    
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      COMMON/JCOM/JBUF(8),JDAT(18528)                                                                                               
      COMMON/ A  / PHI(1871328)                                                                                                     
      COMMON/ B  /   T(1871328)                                                                                                     
                                                                                                                                    
      DATA GZUP /-6.8E04/                                                                                                           
      DATA RGAS /   287.04/                                                                                                         
      DATA MAXX  /  18528/,  MAXL /1871328/,  MAXLEV /100/                                                                          
C---------------------------------------------------------------------                                                              
      NFIL=7                                                                                                                        
      CALL JCLPNT(NFIL, 11,12,13, 14,15, 5,6)                                                                                       
      DO 110 N=11,15                                                                                                                
  110 REWIND N                                                                                                                      
                                                                                                                                    
C     * READ-IN DIRECTIVE CARDS.                                                                                                    
                                                                                                                                    
      READ(5,5010,END=912)  NPL,RLUP,GAMMA,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                                                                                                                         
      WRITE(6,6010) RLUP,GAMMA,ICOORD,PTOIT                                                                                         
      IF(NPL.GT.MAXLEV) CALL                       XIT('GSAPZL ',-1)                                                                
      READ(5,5020,END=913) (LP(I),I=1,NPL)                                      J4                                                  
      CALL LVDCODE(PR,LP,NPL)                                                                                                       
      CALL WRITLEV(PR,NPL,' PR ')                                                                                                   
                                                                                                                                    
      DO 120 L=2,NPL                                                                                                                
  120 IF(PR(L).LE.PR(L-1)) CALL                    XIT('GSAPZL ',-2)                                                                
                                                                                                                                    
      DO 125 L=1,NPL                                                                                                                
  125 PRLOG(L)=LOG(PR(L))                                                                                                           
                                                                                                                                    
C     * GET ETA VALUES FROM PHI FILE.                                                                                               
                                                                                                                                    
      CALL FILEV (LG,NSL,IBUF,11)                                                                                                   
      IF(NSL.LT.1 .OR. NSL.GT.MAXLEV) CALL         XIT('GSAPZL ',-3)                                                                
      CALL LVDCODE(SG,LG,NSL)                                                                                                       
      DO 130 L=1,NSL                                                                                                                
         SG(L) = SG(L)*0.001E0                                                                                                      
  130 CONTINUE                                                                                                                      
      SG(NSL+1) = 1.E0                                                                                                              
      CALL WRITLEV(SG,NSL,' ETA')                                                                                                   
                                                                                                                                    
C     * CONVERT FULL TO HALF ETA LEVELS CONSISTENTLY WITH PAKGCMH.                                                                  
                                                                                                                                    
      DO 140 L=1,NSL                                                                                                                
         SH   (L)=SQRT(SG(L)*SG(L+1))                                                                                               
  140 CONTINUE                                                                                                                      
                                                                                                                                    
C     * DEFINE PARAMETERS OF THE HYBRID VERTICAL DISCRETIZATION.                                                                    
                                                                                                                                    
      CALL COORDAB (AH,BH, NSL  ,SH,ICOORD,PTOIT)                                                                                   
      CALL COORDAB (AG,BG, NSL+1,SG,ICOORD,PTOIT)                                                                                   
                                                                                                                                    
C     * READ MOUNTAINS AND KEEP IN JCOM.                                                                                            
                                                                                                                                    
      CALL RECGET (12,NC4TO8("GRID"),-1,NC4TO8("PHIS"),1,JBUF,MAXX,OK)                                                              
      IF(.NOT.OK) CALL                             XIT('GSAPZL ',-4)                                                                
      WRITE(6,6025) JBUF                                                                                                            
      CALL CMPLBL (0,IBUF, 0,JBUF, OK)                                                                                              
      IF(.NOT.OK) THEN                                                                                                              
         WRITE(6,6025) IBUF,JBUF                                                                                                    
         CALL                                      XIT('GSAPZL ',-5)                                                                
      ENDIF                                                                                                                         
                                                                                                                                    
C     * STOP IF THERE IS NOT ENOUGH SPACE.                                                                                          
                                                                                                                                    
      NWDS = JBUF(5)*JBUF(6)                                                                                                        
      ISP  = MAX(NSL,NPL) *NWDS+1                                                                                                   
      IF((NPL+1)*NWDS.GT.MAXL) CALL                XIT('GSAPZL ',-6)                                                                
      IF((NSL+1)*NWDS.GT.MAXL) CALL                XIT('GSAPZL ',-7)                                                                
C---------------------------------------------------------------------                                                              
C     * GET NEXT SET OF GEOPOTENTIALS INTO ARRAY PHI.                                                                               
                                                                                                                                    
      NSETS=0                                                                                                                       
  150 CALL GETSET2 (11,PHI,LG,ISL,IBUF,MAXX,OK)                                                                                     
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6030) NSETS                                                                                                         
        CALL                                       XIT('GSAPZL ',0)                                                                 
      ENDIF                                                                                                                         
      IF(IBUF(3).NE.NC4TO8(" PHI")) CALL           XIT('GSAPZL ',-8)                                                                
      IF(ISL.NE.NSL) CALL                          XIT('GSAPZL ',-9)                                                                
      CALL CMPLBL(0,IBUF,0,JBUF,OK)                                                                                                 
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6025) IBUF,JBUF                                                                                                     
        CALL                                       XIT('GSAPZL ',-10)                                                               
      ENDIF                                                                                                                         
      NPACK=IBUF(8)                                                                                                                 
                                                                                                                                    
C     * TRANSFER MOUNTAINS FROM JBUF TO LEVEL NSL+1 OF PHI.                                                                         
                                                                                                                                    
      IPHIS=NSL*NWDS+1                                                                                                              
      CALL RECUP2 (PHI(IPHIS),JBUF)                                                                                                 
                                                                                                                                    
C     * GET LN(SF PRES) FOR THIS STEP INTO LAST LEVEL OF T,                                                                         
C     * WHERE THE LAST LEVEL IS THE LARGER OF NSL+1 AND NPL+1.                                                                      
                                                                                                                                    
      NST = IBUF(2)                                                                                                                 
      CALL GETFLD2 (13,T(ISP),NC4TO8("GRID"),NST,NC4TO8("LNSP"),1,                                                                  
     +                                               IBUF,MAXX,OK)                                                                  
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      IF(.NOT.OK) CALL                             XIT('GSAPZL ',-11)                                                               
      CALL CMPLBL(0,IBUF,0,JBUF,OK)                                                                                                 
      IF(.NOT.OK) CALL                             XIT('GSAPZL ',-12)                                                               
                                                                                                                                    
C     * COMPUTE T ON HALF LEVELS FROM "PHI" ON FULL LEVELS IN A WAY THAT                                                            
C     * RECOVERS THE CORRECT TEMPERATURES ON ETA COORDINATES.                                                                       
C     * THEN RECMPUTE THE CORRECT PHI IN ETA COORDINATE,                                                                            
C     * INTERPOLATE T ON ETA LEVELS TO PRESSURE LEVELS, IN-PLACE.                                                                   
                                                                                                                                    
      CALL STFGZ  (T, PHI,NWDS,NSL+1,SG,RGAS)                                                                                       
      CALL ETAPHI (PHI, T,T(ISP),                                                                                                   
     1             AG,BG,RGAS,NWDS,NSL,NSL+1,SIG)                                                                                   
      CALL EAPL   (T, NWDS,PRLOG,NPL,  T,SIG,NSL, T(ISP),RLUP,0.E0,                                                                 
     1             AH,BH,NSL+1,FSIG,DFLNSIG,DLNSIG)                                                                                 
                                                                                                                                    
C     * INTERPOLATE PHI FROM ETA TO PRESSURE, IN-PLACE.                                                                             
C     * EXTRAPOLATION USES ASSUMED DERIVATIVES GZUP,GZDN.                                                                           
C     * THEN CORRECT BOTH T AND GZ UNDER THE GROUND.                                                                                
                                                                                                                                    
      CALL EAPL   (PHI, NWDS,PRLOG,NPL, PHI,SIG,NSL+1, T(ISP),GZUP,0.E0,                                                            
     2             AG,BG,NSL+2,FSIG,DFLNSIG,DLNSIG)                                                                                 
      CALL TGZEX  (PHI, T,T(ISP),NWDS,PR ,NPL,GAMMA)                                                                                
                                                                                                                                    
C     * WRITE-OUT PRESSURE LEVELS PHI AND TEMP.                                                                                     
                                                                                                                                    
      IBUF(3)=NC4TO8(" PHI")                                                                                                        
      IBUF(8)=NPACK                                                                                                                 
      CALL PUTSET2 (14,PHI,LP,NPL,IBUF,MAXX)                                                                                        
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      IBUF(3)=NC4TO8("TEMP")                                                                                                        
      CALL PUTSET2 (15,  T,LP,NPL,IBUF,MAXX)                                                                                        
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      NSETS=NSETS+1                                                                                                                 
      GO TO 150                                                                                                                     
                                                                                                                                    
C     * E.O.F. ON INPUT.                                                                                                            
                                                                                                                                    
  912 CALL                                         XIT('GSAPZL ',-13)                                                               
  913 CALL                                         XIT('GSAPZL ',-14)                                                               
C---------------------------------------------------------------------                                                              
 5010 FORMAT(10X,I5,2E10.0,1X,A4,E10.0)                                         J4                                                  
 5020 FORMAT(16I5)                                                              J4                                                  
 6010 FORMAT('0 RLUP,GAMMA =',2E12.3,', ICOORD=',1X,A4,                                                                             
     1       ', P.LID (PA)=',E10.3)                                                                                                 
 6025 FORMAT(' ',A4,I10,2X,A4,I10,4I6)                                                                                              
 6030 FORMAT('0',I5,' SETS INTERPOLATED')                                                                                           
      END                                                                                                                           
