      PROGRAM TAPHI                                                                                                                 
C     PROGRAM TAPHI (GSTEMP,       GSPHIS,       GSLNSP,       GSRGASM,         J2                                                  
C    1                             GSPHI,         INPUT,        OUTPUT, )       J2                                                  
C    2        TAPE11=GSTEMP,TAPE12=GSPHIS,TAPE13=GSLNSP,TAPE14=GSRGASM,                                                             
C    3                      TAPE15=GSPHI , TAPE5= INPUT,TAPE6 = OUTPUT)                                                             
C     -----------------------------------------------------------------         J2                                                  
C                                                                               J2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       J2                                                  
C     JUL 27/94 - J. KOSHYK                                                                                                         
C                                                                               J2                                                  
CTAPHI   - CALCULATES PHI ON ETA (SIGMA/HYBRID) LEVELS FROM TEMPERATURE         J1                                                  
C          AND MOIST GAS CONSTANT ON ETA LEVELS.                        4  1 C  J1                                                  
C                                                                               J3                                                  
CAUTHOR  - J. KOSHYK                                                            J3                                                  
C                                                                               J3                                                  
CPURPOSE - CALCULATES PHI ON ETA (SIGMA/HYBRID) LEVELS FROM TEMPERATURE         J3                                                  
C          AND MOIST GAS CONSTANT ON ETA LEVELS.                                J3                                                  
CINPUT FILES...                                                                 J3                                                  
C                                                                               J3                                                  
C      GSTEMP  = SETS OF TEMPERATURE ON ETA LEVELS.                             J3                                                  
C      GSPHIS  = MOUNTAINS (SURFACE GEOPOTENTIAL).                              J3                                                  
C      GSLNSP  = SETS OF LN(SURFACE PRESSURE).                                  J3                                                  
C      GSRGASM = SETS OF MOIST GAS CONSTANT ON ETA LEVELS.                      J3                                                  
C                                                                               J3                                                  
COUTPUT FILES...                                                                J3                                                  
C                                                                               J3                                                  
C      GSPHI   = ETA LEVEL GEOPOTENTIALS (ON SAME LEVELS AS GSTEMP).            J3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               J5                                                  
C      ICOORD = 4H SIG/4H ETA/4HET10/4HET15  FOR VERTICAL COORDINATE.           J5                                                  
C      LAY    = LAYERING SCHEME FOR CALCULATING LAYER INTERFACES.               J5                                                  
C      PTOIT  = PRESSURE (PA) OF LID AT TOP OF MODEL.                           J5                                                  
C                                                                               J5                                                  
CEXAMPLE OF INPUT CARDS...                                                      J5                                                  
C                                                                               J5                                                  
C*TAPHI.    ET15    3        0.                                                 J5                                                  
C-------------------------------------------------------------------------                                                          
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
                                                                                                                                    
      INTEGER LH(101)                                                                                                               
                                                                                                                                    
      REAL SFM(101),SHM(100),SFT(101),SHT(100)                                                                                      
      REAL AFT(101),BFT(101),AHT(100),BHT(100)                                                                                      
      REAL SIGFT(101),SIGHT(100),DLNSIGFT(100)                                                                                      
      REAL RGASM(1871328)                                                                                                           
                                                                                                                                    
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      COMMON/JCOM/JBUF(8),JDAT(18528)                                                                                               
      COMMON/ A  / PHI(1871328)                                                                                                     
      COMMON/ B  /   T(1871328)                                                                                                     
                                                                                                                                    
      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) ICOORD,LAY,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) ICOORD,LAY,PTOIT                                                                                                
                                                                                                                                    
C     * GET ETA VALUES FROM TEMPERATURE FILE.                                                                                       
C     * INITIALIZE SHM ARRAY TO SHT IN ORDER TO USE BASCAL.                                                                         
C     * (NOTE: SFM = FULL MOMENTUM LAYERS,                                                                                          
C     *        SHM = HALF MOMENTUM LAYERS,                                                                                          
C     *        SFT = FULL THERMOD. LAYERS,                                                                                          
C     *        SHT = HALF THERMOD. LAYERS.)                                                                                         
                                                                                                                                    
      CALL FILEV (LH,NSL,IBUF,11)                                                                                                   
      IF(NSL.LT.1 .OR. NSL.GT.MAXLEV) CALL         XIT('TAPHI ',-1)                                                                 
      CALL LVDCODE(SHT,LH,NSL)                                                                                                      
      DO 130 L=1,NSL                                                                                                                
         SHT(L) = SHT(L)*0.001E0                                                                                                    
         SHM(L) = SHT(L)                                                                                                            
  130 CONTINUE                                                                                                                      
      CALL WRITLEV(SHT,NSL,' ETA')                                                                                                  
                                                                                                                                    
C     * COMPUTE INTERFACE ("FULL") LEVELS FOR PHI BASED ON LAYERING                                                                 
C     * PARAMETER.                                                                                                                  
                                                                                                                                    
      CALL BASCAL(SFM,SFT(2),SHM,SHT,NSL,LAY)                                                                                       
      SFT(1)=PTOIT/101320.E0                                                                                                        
                                                                                                                                    
C     * CALCULATE PARAMETERS OF THE HYBRID VERTICAL DISCRETIZATION,                                                                 
C     * A AND B, FOR ALL THERMOD. LAYERS (HALF AND FULL).                                                                           
                                                                                                                                    
      CALL COORDAB (AHT,BHT, NSL  ,SHT ,ICOORD,PTOIT)                                                                               
      CALL COORDAB (AFT,BFT, NSL+1,SFT ,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('TAPHI ',-2)                                                                 
      WRITE(6,6025) JBUF                                                                                                            
      CALL CMPLBL (0,IBUF, 0,JBUF, OK)                                                                                              
      IF(.NOT.OK) THEN                                                                                                              
         WRITE(6,6025) IBUF,JBUF                                                                                                    
         CALL                                      XIT('TAPHI ',-3)                                                                 
      ENDIF                                                                                                                         
                                                                                                                                    
C     * STOP IF THERE IS NOT ENOUGH SPACE.                                                                                          
                                                                                                                                    
      NWDS = JBUF(5)*JBUF(6)                                                                                                        
      IF((NSL+1)*NWDS.GT.MAXL) CALL                XIT('TAPHI ',-4)                                                                 
C---------------------------------------------------------------------                                                              
C     * GET NEXT SET OF TEMPERATURES INTO ARRAY T.                                                                                  
                                                                                                                                    
      NSETS=0                                                                                                                       
  150 CALL GETSET2 (11,T,LH,ISL,IBUF,MAXX,OK)                                                                                       
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6030) NSETS                                                                                                         
        CALL                                       XIT('TAPHI  ',0)                                                                 
      ENDIF                                                                                                                         
      IF(IBUF(3).NE.NC4TO8("TEMP")) CALL           XIT('TAPHI ',-5)                                                                 
      IF(ISL.NE.NSL) CALL                          XIT('TAPHI ',-6)                                                                 
      CALL CMPLBL(0,IBUF,0,JBUF,OK)                                                                                                 
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6025) IBUF,JBUF                                                                                                     
        CALL                                       XIT('TAPHI ',-7)                                                                 
      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     * PUT LN(SF PRES) FOR THIS STEP INTO LEVEL NSL+1 OF T.                                                                        
                                                                                                                                    
      ILNSP=IPHIS                                                                                                                   
      NST = IBUF(2)                                                                                                                 
      CALL GETFLD2 (13,T(ILNSP),NC4TO8("GRID"),NST,NC4TO8("LNSP"),                                                                  
     +                                             1,IBUF,MAXX,OK)                                                                  
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      IF(.NOT.OK) CALL                             XIT('TAPHI ',-8)                                                                 
      CALL CMPLBL(0,IBUF,0,JBUF,OK)                                                                                                 
      IF(.NOT.OK) CALL                             XIT('TAPHI ',-9)                                                                 
                                                                                                                                    
C     * GET NEXT SET OF MOIST GAS CONSTANT INTO ARRAY RGASM.                                                                        
                                                                                                                                    
      CALL GETSET2 (14,RGASM,LH,ISL,IBUF,MAXX,OK)                                                                                   
      IF(NSETS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      IF(.NOT.OK) CALL                             XIT('TAPHI ',-10)                                                                
      IF(IBUF(3).NE.NC4TO8("RGAS")) CALL           XIT('TAPHI ',-11)                                                                
      IF(ISL.NE.NSL) CALL                          XIT('TAPHI ',-12)                                                                
      CALL CMPLBL(0,IBUF,0,JBUF,OK)                                                                                                 
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6025) IBUF,JBUF                                                                                                     
        CALL                                       XIT('TAPHI ',-13)                                                                
      ENDIF                                                                                                                         
                                                                                                                                    
C     * COMPUTE PHI ON MODEL HALF LEVELS USING HYDROSTATIC EQUATION.                                                                
                                                                                                                                    
      CALL PHIETA (PHI, T,T(ILNSP), AFT,BFT,AHT,BHT,                                                                                
     1             RGASM, NWDS,NSL,NSL+1, SIGFT,SIGHT,DLNSIGFT)                                                                     
                                                                                                                                    
C     * WRITE-OUT PHI ON ETA LEVELS.                                                                                                
                                                                                                                                    
      IBUF(3)=NC4TO8(" PHI")                                                                                                        
      IBUF(8)=NPACK                                                                                                                 
      CALL PUTSET2 (15,PHI,LH,NSL,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('TAPHI ',-14)                                                                
  913 CALL                                         XIT('TAPHI ',-15)                                                                
C---------------------------------------------------------------------------------                                                  
 5010 FORMAT(10X,1X,A4,I5,E10.0)                                                J4                                                  
 6010 FORMAT('0 ICOORD=',1X,A4,', LAY=',I5,', P.LID (PA)=',E10.3)                                                                   
 6025 FORMAT(' ',A4,I10,2X,A4,1X,I10,4I6)                                                                                           
 6030 FORMAT('0',I5,' SETS INTERPOLATED')                                                                                           
      END                                                                                                                           
