      PROGRAM GSTMTD                                                                                                                
C     PROGRAM GSTMTD (GSES,       GSTEMP,       GSLNSP,       GSTSTD,           I2                                                  
C    1                                          INPUT,        OUTPUT,   )       I2                                                  
C    2          TAPE1=GSES, TAPE2=GSTEMP, TAPE3=GSLNSP, TAPE8=GSTSTD,                                                               
C    3                                    TAPE5=INPUT,  TAPE6=OUTPUT)                                                               
C     ---------------------------------------------------------------           I2                                                  
C                                                                               I2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       I2                                                  
C     JAN 13/93 - E. CHAN     (DECODE LEVELS IN 8-WORD LABEL)                                                                       
C     JAN 29/92 - E. CHAN     (CONVERT HOLLERITH LITERALS TO ASCII)                                                                 
C     DEC 18/90 - M.LAZARE. - ADD 4HQHYB AS ANOTHER CHOICE FOR MOIST.                                                               
C     MAR 23/89 - M.LAZARE.                                                                                                         
C                                                                               I2                                                  
CGSTMTD  - CONVERTS ETA (SIGMA/HYBRID) LEVEL MOISTURE VARIABLE TO               I1                                                  
C          DEWPOINT DEPRESSION FOR MODEL INITIALIZATION.                3  2 C  I1                                                  
C                                                                               I3                                                  
CAUTHOR  - M.LAZARE                                                             I3                                                  
C                                                                               I3                                                  
CPURPOSE - CONVERTS ETA (SIGMA/HYBRID) LEVEL GRID FILE OF MODEL                 I3                                                  
C          MOISTURE VARIABLE TO DEWPOINT DEPRESSION FOR MODEL                   I3                                                  
C          INITIALIZATION.                                                      I3                                                  
C                                                                               I3                                                  
CINPUT FILES...                                                                 I3                                                  
C                                                                               I3                                                  
C      GSES   = MODEL MOISTURE VARIABLE (ES) ON ETA (SIGMA/HYBRID) LEVELS.      I3                                                  
C      GSTEMP = TEMPERATURES ON THE SAME (OR MORE) LEVELS.                      I3                                                  
C      GSLNSP = CORRESPONDING SERIES OF LN(SF PRES), IN MB.                     I3                                                  
C                                                                               I3                                                  
COUTPUT FILES...                                                                I3                                                  
C                                                                               I3                                                  
C      GSTSTD = OUTPUT FILE FOR THE DEWPOINT DEPRESSION.                        I3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               I5                                                  
C      ICOORD = 4H SIG/4H ETA FOR SIGMA/ETA VERTICAL COORDINATE.                I5                                                  
C      MOIST  = 4HT-TD/4H  TD/4HRLNQ/4HSQRT/4H   Q FOR VARIABLE                 I5                                                  
C               DEW POINT DEPRESSION/DEW POINT/-(LN(Q)**-1)/Q**0.5/Q            I5                                                  
C      PTOIT  = PRESSURE (PA) OF THE RIGID LID OF THE MODEL.                    I5                                                  
C                                                                               I5                                                  
CEXAMPLE OF INPUT CARD...                                                       I5                                                  
C                                                                               I5                                                  
C*GSTMTD.   SIG    Q        0.                                                  I5                                                  
C------------------------------------------------------------------------                                                           
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
                                                                                                                                    
      REAL ETAH(100),AH(100),BH(100)                                                                                                
                                                                                                                                    
      INTEGER LETAT(100),LETAS(100),KBUF(8)                                                                                         
                                                                                                                                    
      COMMON /LEVELS/ SIGH  (18528)                                                                                                 
      COMMON /BLANCK/ PSMBLN(18528),F   (18528),G   (18528)                                                                         
      REAL            PRESS (18528)                                                                                                 
      EQUIVALENCE    (PRESS,PSMBLN)                                                                                                 
                                                                                                                                    
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      COMMON/EPS / A, B, EPS1, EPS2                                                                                                 
                                                                                                                                    
      DATA MAXX/18528/, MAXLEV/100/                                                                                                 
C---------------------------------------------------------------------                                                              
      NFIL = 6                                                                                                                      
      CALL JCLPNT (NFIL,1,2,3,8,5,6)                                                                                                
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      REWIND 3                                                                                                                      
      REWIND 8                                                                                                                      
C                                                                                                                                   
C     * GET THERMODYNAMIC CONSTANTS.                                                                                                
C                                                                                                                                   
      CALL SPWCONH(FVORT,PI)                                                                                                        
                                                                                                                                    
C     * READ-IN DIRECTIVE CARD.                                                                                                     
                                                                                                                                    
      READ(5,5000,END=911) ICOORD,MOIST,PTOIT                                   I4                                                  
      IF(ICOORD.EQ.NC4TO8("    ")) ICOORD=NC4TO8(" SIG")                                                                            
      IF( MOIST.EQ.NC4TO8("    "))  MOIST=NC4TO8("T-TD")                                                                            
      IF(ICOORD.EQ.NC4TO8(" SIG")) THEN                                                                                             
        PTOIT=MAX(PTOIT,0.00E0)                                                                                                     
      ELSE                                                                                                                          
        PTOIT=MAX(PTOIT, 0.100E-09)                                                                                                 
      ENDIF                                                                                                                         
      WRITE(6,6000) ICOORD,MOIST,PTOIT                                                                                              
                                                                                                                                    
C     * GET ETA VALUES FROM TEMP FILE. ALSO FROM ES FILE.                                                                           
C     * ILEV,LEVS = NUMBER OF TEMP,MOISTURE LEVELS IN THE MODEL.                                                                    
C     * ETAH = ETA VALUES OF THE TEMPERATURE LEVELS.                                                                                
C     * MISSING MOISTURE LEVELS ARE AT THE TOP (LEVELS 1,2,..ETC).                                                                  
                                                                                                                                    
      CALL FILEV (LETAT ,ILEV,IBUF,2)                                                                                               
      IF(ILEV.LT.1 .OR. ILEV.GT.MAXLEV) CALL       XIT('GSTMTD',-1)                                                                 
      WRITE(6,6015) IBUF(3),ILEV,(LETAT(L),L=1,ILEV)                                                                                
      DO 116 I=1,8                                                                                                                  
  116 KBUF(I)=IBUF(I)                                                                                                               
      CALL LVDCODE(ETAH,LETAT,ILEV)                                                                                                 
      DO 120 L=1,ILEV                                                                                                               
  120 ETAH(L)=ETAH(L)*0.001E0                                                                                                       
                                                                                                                                    
      CALL FILEV (LETAS,LEVS,IBUF,1)                                                                                                
      IF(LEVS.LT.1 .OR. LEVS.GT.ILEV) CALL         XIT('GSTMTD',-2)                                                                 
      WRITE(6,6015) IBUF(3),LEVS,(LETAS(L),L=1,LEVS)                                                                                
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSTMTD',-3)                                                                 
      MISSING=ILEV-LEVS                                                                                                             
      NWDS   =IBUF(5)*IBUF(6)                                                                                                       
                                                                                                                                    
C     * DEFINE PARAMETERS OF THE VERTICAL DISCRETIZATION.                                                                           
                                                                                                                                    
      CALL COORDAB (AH,BH, ILEV,ETAH,ICOORD,PTOIT)                                                                                  
C---------------------------------------------------------------------                                                              
C     * GET THE NEXT LNSP AND CONVERT TO SFC.PRES. IN PA.                                                                           
                                                                                                                                    
      NRECS=0                                                                                                                       
  200 CALL GETFLD2 (3,PSMBLN,NC4TO8("GRID"),-1,NC4TO8("LNSP"),1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF(NRECS.EQ.0) WRITE(6,6025) IBUF                                                                                             
      IF(.NOT.OK)THEN                                                                                                               
        WRITE(6,6010) NRECS                                                                                                         
        IF(NRECS.EQ.0)THEN                                                                                                          
          CALL                                     XIT('GSTMTD',-4)                                                                 
        ELSE                                                                                                                        
          CALL                                     XIT('GSTMTD',0)                                                                  
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      ITIM=IBUF(2)                                                                                                                  
                                                                                                                                    
      DO 210 I=1,NWDS                                                                                                               
         PRESS(I)=100.E0*EXP(PSMBLN(I))                                                                                             
  210 CONTINUE                                                                                                                      
                                                                                                                                    
C     * LEVEL LOOP OVER ALL TEMPERATURE LEVELS.                                                                                     
C       --------------------------------------                                                                                      
                                                                                                                                    
      DO 500 L= 1,ILEV                                                                                                              
                                                                                                                                    
C        * DEFINE LOCAL SIGMA VALUES.                                                                                               
                                                                                                                                    
         CALL NIVCAL (SIGH, AH(L),BH(L),PRESS,1,NWDS,NWDS)                                                                          
                                                                                                                                    
C        * GET THE TEMPERATURE INTO ARRAY G.                                                                                        
                                                                                                                                    
         CALL GETFLD2 (2,G,NC4TO8("GRID"),ITIM,NC4TO8("TEMP"),                                                                      
     +                                  LETAT(L),IBUF,MAXX,OK)                                                                      
         IF(NRECS.EQ.0) WRITE(6,6025) IBUF                                                                                          
         IF(.NOT.OK)THEN                                                                                                            
           WRITE(6,6010) NRECS                                                                                                      
           CALL                                    XIT('GSTMTD',-5)                                                                 
         ENDIF                                                                                                                      
         CALL CMPLBL (0,IBUF, 0,KBUF, OK)                                                                                           
         IF(.NOT.OK) CALL                          XIT('GSTMTD',-6)                                                                 
                                                                                                                                    
C        * IF THERE ARE MISSING MOISTURE LEVELS AT THE TOP                                                                          
C        * JUST SET T-TD TO 100.                                                                                                    
                                                                                                                                    
         IF(L.LE.MISSING) THEN                                                                                                      
            DO 220 I=1,NWDS                                                                                                         
              F(I)=100.E0                                                                                                           
  220       CONTINUE                                                                                                                
         ELSE                                                                                                                       
                                                                                                                                    
C           * GET THE MODEL MOISTURE VARIABLE INTO ARRAY F, THEN                                                                    
C           * CALCULATE T-TD.                                                                                                       
                                                                                                                                    
            CALL GETFLD2 (1,F,NC4TO8("GRID"),ITIM,NC4TO8("  ES"),                                                                   
     +                                     LETAT(L),IBUF,MAXX,OK)                                                                   
            IF(NRECS.EQ.0) WRITE(6,6025) IBUF                                                                                       
            IF(.NOT.OK)THEN                                                                                                         
              WRITE(6,6010) NRECS                                                                                                   
              CALL                                 XIT('GSTMTD',-7)                                                                 
            ENDIF                                                                                                                   
            CALL CMPLBL (0,IBUF, 0,KBUF, OK)                                                                                        
            IF(.NOT.OK) CALL                       XIT('GSTMTD',-8)                                                                 
                                                                                                                                    
            IF(MOIST.EQ.NC4TO8("T-TD")) THEN                                                                                        
               CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("  TD")) THEN                                                                                    
               DO 315 I=1,NWDS                                                                                                      
                  F(I )= G(I) - F(I)                                                                                                
  315          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("RLNQ"))THEN                                                                                     
               DO 320 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  ES   = F(I)                                                                                                       
                  SHUM = EXP(-1.E0/ES)                                                                                              
                  ALVP = LOG (SHUM*PRES/(EPS1+EPS2*SHUM) )                                                                          
                  TD   = B/(A - ALVP)                                                                                               
                  F(I) = G(I) - TD                                                                                                  
  320          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("   Q"))THEN                                                                                     
               DO 325 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  SHUM = MAX(F(I),1.E-20)                                                                                           
                  ALVP = LOG (SHUM*PRES/(EPS1+EPS2*SHUM) )                                                                          
                  TD   = B/(A - ALVP)                                                                                               
                  F(I) = G(I) - TD                                                                                                  
  325          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("SQRT"))THEN                                                                                     
               DO 330 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  ES   = F(I)                                                                                                       
                  SHUM = MAX(ES**2,1.E-20)                                                                                          
                  ALVP = LOG (SHUM*PRES/(EPS1+EPS2*SHUM) )                                                                          
                  TD   = B/(A - ALVP)                                                                                               
                  F(I) = G(I) - TD                                                                                                  
  330          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8(" LNQ"))THEN                                                                                     
               DO 335 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  ES   = F(I)                                                                                                       
                  SHUM = EXP(ES)                                                                                                    
                  ALVP = LOG (SHUM*PRES/(EPS1+EPS2*SHUM) )                                                                          
                  TD   = B/(A - ALVP)                                                                                               
                  F(I) = G(I) - TD                                                                                                  
  335          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("QHYB"))THEN                                                                                     
               SHUMREF=10.E-3                                                                                                       
               QMIN=1.E-20                                                                                                          
               SMIN=SHUMREF/(1.E0+LOG(SHUMREF/QMIN))                                                                                
               DO 340 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  ES   = F(I)                                                                                                       
                  IF(ES.GE.SHUMREF) THEN                                                                                            
                     SHUM=ES                                                                                                        
                  ELSE                                                                                                              
                     ES = MAX(ES,SMIN)                                                                                              
                     SHUM   = SHUMREF/(EXP(SHUMREF/ES-1.E0))                                                                        
                  ENDIF                                                                                                             
                  ALVP = LOG (SHUM*PRES/(EPS1+EPS2*SHUM) )                                                                          
                  TD   = B/(A - ALVP)                                                                                               
                  F(I) = G(I) - TD                                                                                                  
  340          CONTINUE                                                                                                             
                                                                                                                                    
            ELSE                                                                                                                    
               CALL                                XIT('GSTMTD',-9)                                                                 
            ENDIF                                                                                                                   
         ENDIF                                                                                                                      
                                                                                                                                    
C        * BOUND T-TD TO ENSURE INITIALIZATION WILL WORK PROPERLY.                                                                  
                                                                                                                                    
         DO 350 I=1,NWDS                                                                                                            
            F(I) = MAX(0.E0, MIN(F(I),100.E0) )                                                                                     
  350    CONTINUE                                                                                                                   
                                                                                                                                    
C        * PUT T-TD ON FILE 8.                                                                                                      
                                                                                                                                    
         IBUF(4)=LETAT(L)                                                                                                           
         IBUF(3)=NC4TO8("  ES")                                                                                                     
         CALL PUTFLD2 (8,F,IBUF,MAXX)                                                                                               
         IF(NRECS.EQ.0) WRITE(6,6025) IBUF                                                                                          
                                                                                                                                    
         NRECS=NRECS+1                                                                                                              
                                                                                                                                    
  500 CONTINUE                                                                                                                      
                                                                                                                                    
      GO TO 200                                                                                                                     
                                                                                                                                    
C     * E.O.F. ON INPUT.                                                                                                            
  911 CALL                                         XIT('GSTMTD',-10)                                                                
C---------------------------------------------------------------------                                                              
 5000 FORMAT(10X,2(1X,A4),E10.0)                                                I4                                                  
 6000 FORMAT(' COORD = ',A4,', MOISTURE VARIABLE = ',A4,                                                                            
     1       ', P.LID (PA)=',E10.3)                                                                                                 
 6010 FORMAT('0 GSTMTD TRANSFORMED',I6, '  PAIRS')                                                                                  
 6015 FORMAT(' NAME=',A4,',ILEV=',I5/(' ETA LEVELS=',10I6)/(12X,10I6))                                                              
 6025 FORMAT(' ',A4,I10,2X,A4,5I10)                                                                                                 
      END                                                                                                                           
