      PROGRAM GSHUMH                                                                                                                
C     PROGRAM GSHUMH (GSES,       GSTEMP,       GSLNSP,       GSSHUM,           J2                                                  
C    1                            GSRHUM,       INPUT,        OUTPUT,   )       J2                                                  
C    2          TAPE1=GSES, TAPE2=GSTEMP, TAPE3=GSLNSP, TAPE8=GSSHUM,                                                               
C    3                      TAPE9=GSRHUM, TAPE5=INPUT,  TAPE6=OUTPUT)                                                               
C     ---------------------------------------------------------------           J2                                                  
C                                                                               J2                                                  
C     DEC 07/09 - S.KHARIN (REVISED QMIN=MAX(SREF*1.E-16,1.E-36)                J2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                           
C     JUN 10/02 - M.LAZARE. - GENERALIZED READ-IN OF SREF (FORMERLY                                                                 
C                             SHUMREF) AND SPOW FOR QHYB,SLQB.                                                                      
C     APR 10/00 - M.LAZARE. - ADD 4HSLQB AS ANOTHER CHOICE FOR MOIST.                                                               
C     JUL 22/96 - M.LAZARE. - ADD 4HSL3D AS ANOTHER CHOICE FOR MOIST.                                                               
C     DEC 22/94 - F.MAJAESS (ALLOW READING THE REFERENCE VALUE)                                                                     
C     JAN 12/93 - E. CHAN   (DECODE LEVELS IN 8-WORD LABEL)                                                                         
C     JAN 29/92 - E. CHAN   (CONVERT HOLLERITH LITERALS TO ASCII)                                                                   
C     MAY 17/90 - M.LAZARE. - ADD 4HQHYB AS ANOTHER CHOICE FOR MOIST.                                                               
C     JUN 15/89 - M.LAZARE. - ADD 4H LNQ AS ANOTHER CHOICE FOR MOIST.                                                               
C     AUG 11/88 - M.LAZARE  (REPLACE "A,B,EPS1,EPS2" DATA BY COMMON STATEMENT.                                                      
C                            ALSO, INSERT A CALL TO SPWCONH ROUTINE)                                                                
C     FEB 29/88 - R.LAPRISE.                                                                                                        
C                                                                               J2                                                  
CGSHUMH  - CONVERTS ETA (SIGMA/HYBRID) LEVEL MOISTURE VARIABLE TO               J1                                                  
C          SPECIFIC AND RELATIVE HUMIDITY                               3  2 C  J1                                                  
C                                                                               J3                                                  
CAUTHOR - R. LAPRISE                                                            J3                                                  
C                                                                               J3                                                  
CPURPOSE - CONVERTS ETA (SIGMA/HYBRID) LEVEL GRID FILE OF MODEL                 J3                                                  
C          MOISTURE VARIABLE TO SPECIFIC AND RELATIVE HUMIDITY.                 J3                                                  
C          NOTE - THIS PROGRAM IS UPWARD COMPATIBLE WITH THE                    J3                                                  
C                 EARLIER VERSION (GSHUM PROGRAM), EXCEPT THAT                  J3                                                  
C                 IT NOW READS AN INPUT CARD.                                   J3                                                  
C                                                                               J3                                                  
CINPUT FILES...                                                                 J3                                                  
C                                                                               J3                                                  
C      GSES   = MODEL MOISTURE VARIABLE (ES) ON ETA (SIGMA/HYBRID) LEVELS.      J3                                                  
C      GSTEMP = TEMPERATURES ON THE SAME (OR MORE) LEVELS.                      J3                                                  
C      GSLNSP = CORRESPONDING SERIES OF LN(SF PRES), IN MB.                     J3                                                  
C                                                                               J3                                                  
COUTPUT FILES...                                                                J3                                                  
C                                                                               J3                                                  
C      GSSHUM = OUTPUT FILE FOR THE SPECIFIC HUMIDITY.                          J3                                                  
C      GSRHUM = OUTPUT FILE FOR THE RELATIVE HUMIDITY.                          J3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               J5                                                  
C      ICOORD = 4H SIG/4H ETA FOR SIGMA/ETA VERTICAL COORDINATE.                J5                                                  
C      MOIST  = 4HT-TD  FOR  VARIABLE DEW POINT DEPRESSION                      J5                                                  
C               4H  TD                DEW POINT                                 J5                                                  
C               4HRLNQ                -(LN(Q)**-1)                              J5                                                  
C               4HSQRT                Q**0.5                                    J5                                                  
C               4H   Q                Q                                         J5                                                  
C               4HSL3D                SL3D (SEMI-LAGRANGIAN OPTION)             J5                                                  
C               4H LNQ                LN(Q)                                     J5                                                  
C               4HQHYB                HYBRID ( =Q                      ; Q<=Q0  J5                                                  
C                                              =Q0/((1+P*LN(Q0/Q))**P  ; Q>Q0   J5                                                  
C                                              WHERE Q0 IS THE REFERENCE VALUE  J5                                                  
C                                              "SREF" BELOW AND "P" IS THE      J5                                                  
C                                              REFERENCE POWER ("SPOW" BELOW)   J5                                                  
C               4HSLQB                QHYB AS ABOVE EXCEPT SEMI-LAGRANGIAN ADV. J5                                                  
C      PTOIT  = PRESSURE (PA) OF THE RIGID LID OF THE MODEL.                    J5                                                  
C      SPOW= REFERENCE POWER FOR THE HYBRID CASE (DEFAULT; 1.)                  J5                                                  
C      SREF= REFERENCE VALUE FOR THE HYBRID CASE (DEFAULT; 10.E-3 KG/KG)        J5                                                  
C                                                                               J5                                                  
CEXAMPLE OF INPUT CARD...                                                       J5                                                  
C                                                                               J5                                                  
C*GSHUMH.   SIG T-TD        0.                                                  J5                                                  
C------------------------------------------------------------------------                                                           
C                                                                                                                                   
C     * SET IN "PTMIN" PLID MINIMUM VALUE ...                                                                                       
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      REAL PTMIN                                                                                                                    
      PARAMETER (PTMIN= 0.100E-09)                                                                                                  
                                                                                                                                    
      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 = 7                                                                                                                      
      CALL JCLPNT (NFIL, 1,2,3, 8,9, 5,6)                                                                                           
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      REWIND 3                                                                                                                      
      REWIND 8                                                                                                                      
      REWIND 9                                                                                                                      
C                                                                                                                                   
C     * GET THERMODYNAMIC CONSTANTS.                                                                                                
C                                                                                                                                   
      CALL SPWCON7(FVORT,PI)                                                                                                        
                                                                                                                                    
C     * READ-IN DIRECTIVE CARD.                                                                                                     
                                                                                                                                    
      READ(5,5000,END=911) ICOORD,MOIST,PTOIT,SREF,SPOW                         J4                                                  
      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,PTMIN)                                                                                                      
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * SET DEFAULT VALUES.                                                                                                         
C                                                                                                                                   
      IF(SREF.LE.0.E0) SREF=10.E-3                                                                                                  
      IF(SPOW.LE.0.E0) SPOW=1.E0                                                                                                    
                                                                                                                                    
      WRITE(6,6000) ICOORD,MOIST,PTOIT                                                                                              
      WRITE(6,6005) SREF,SPOW                                                                                                       
                                                                                                                                    
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.EQ.0 .OR. ILEV.GT.MAXLEV) CALL       XIT('GSHUMH',-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.EQ.0 .OR. LEVS.GT.ILEV) CALL         XIT('GSHUMH',-2)                                                                 
      WRITE(6,6015) IBUF(3),LEVS,(LETAS(L),L=1,LEVS)                                                                                
      CALL CMPLBL (0,IBUF,0,KBUF,OK)                                                                                                
      IF(.NOT.OK) CALL                             XIT('GSHUMH',-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('GSHUMH',-4)                                                                 
        ELSE                                                                                                                        
          CALL                                     XIT('GSHUMH',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('GSHUMH',-5)                                                                 
         ENDIF                                                                                                                      
         CALL CMPLBL (0,IBUF, 0,KBUF, OK)                                                                                           
         IF(.NOT.OK) CALL                          XIT('GSHUMH',-6)                                                                 
                                                                                                                                    
C        * IF THERE ARE MISSING MOISTURE LEVELS AT THE TOP                                                                          
C        * JUST SET RHUM AND SHUM TO ZERO.                                                                                          
                                                                                                                                    
         IF(L.LE.MISSING) THEN                                                                                                      
            DO 220 I=1,NWDS                                                                                                         
              F(I)=0.E0                                                                                                             
              G(I)=0.E0                                                                                                             
  220       CONTINUE                                                                                                                
                                                                                                                                    
C           * GET THE MODEL MOISTURE VARIABLE INTO ARRAY F, THEN                                                                    
C           * CALCULATE SPEC. HUM. IN F, AND REL. HUM. IN G.                                                                        
                                                                                                                                    
         ELSE                                                                                                                       
            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('GSHUMH',-7)                                                                 
            ENDIF                                                                                                                   
            CALL CMPLBL (0,IBUF, 0,KBUF, OK)                                                                                        
            IF(.NOT.OK) CALL                       XIT('GSHUMH',-8)                                                                 
                                                                                                                                    
            IF(MOIST.EQ.NC4TO8("T-TD")) THEN                                                                                        
               DO 310 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  T    = G(I)                                                                                                       
                  ES   = F(I)                                                                                                       
                  TD   = T-ES                                                                                                       
                  E    = EXP(A-B/TD)                                                                                                
                  SHUM = EPS1*E/(PRES-EPS2*E)                                                                                       
                  F(I )= SHUM                                                                                                       
                  E    = EXP(A-B/T)                                                                                                 
                  SAT  = EPS1*E/(PRES-EPS2*E)                                                                                       
                  G(I) = SHUM/SAT                                                                                                   
  310          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("  TD")) THEN                                                                                    
               DO 315 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  T    = G(I)                                                                                                       
                  ES   = F(I)                                                                                                       
                  TD   = ES                                                                                                         
                  E    = EXP(A-B/TD)                                                                                                
                  SHUM = EPS1*E/(PRES-EPS2*E)                                                                                       
                  F(I )= SHUM                                                                                                       
                  E    = EXP(A-B/T)                                                                                                 
                  SAT  = EPS1*E/(PRES-EPS2*E)                                                                                       
                  G(I) = SHUM/SAT                                                                                                   
  315          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("RLNQ"))THEN                                                                                     
               DO 320 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  T    = G(I)                                                                                                       
                  ES   = F(I)                                                                                                       
                  SHUM = EXP(-1.E0/ES)                                                                                              
                  F(I )= SHUM                                                                                                       
                  E    = EXP(A-B/T)                                                                                                 
                  SAT  = EPS1*E/(PRES-EPS2*E)                                                                                       
                  G(I) = SHUM/SAT                                                                                                   
  320          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("   Q") .OR.                                                                                     
     +             MOIST.EQ.NC4TO8("SL3D")     ) THEN                                                                               
               DO 325 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  T    = G(I)                                                                                                       
                  ES   = F(I)                                                                                                       
                  SHUM = ES                                                                                                         
                  F(I )= SHUM                                                                                                       
                  E    = EXP(A-B/T)                                                                                                 
                  SAT  = EPS1*E/(PRES-EPS2*E)                                                                                       
                  G(I) = SHUM/SAT                                                                                                   
  325          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("SQRT"))THEN                                                                                     
               DO 330 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  T    = G(I)                                                                                                       
                  ES   = F(I)                                                                                                       
                  SHUM = ES**2                                                                                                      
                  F(I )= SHUM                                                                                                       
                  E    = EXP(A-B/T)                                                                                                 
                  SAT  = EPS1*E/(PRES-EPS2*E)                                                                                       
                  G(I) = SHUM/SAT                                                                                                   
  330          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8(" LNQ"))THEN                                                                                     
               DO 335 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  T    = G(I)                                                                                                       
                  ES   = F(I)                                                                                                       
                  SHUM = EXP(ES)                                                                                                    
                  F(I )= SHUM                                                                                                       
                  E    = EXP(A-B/T)                                                                                                 
                  SAT  = EPS1*E/(PRES-EPS2*E)                                                                                       
                  G(I) = SHUM/SAT                                                                                                   
  335          CONTINUE                                                                                                             
                                                                                                                                    
            ELSEIF(MOIST.EQ.NC4TO8("QHYB") .OR.                                                                                     
     +             MOIST.EQ.NC4TO8("SLQB")     ) THEN                                                                               
               PINV=1.E0/SPOW                                                                                                       
               QMIN=MAX(SREF*1.E-16,1.E-36)                                                                                         
               SMIN   = SREF/((1.E0+SPOW*LOG(SREF/QMIN))**PINV)                                                                     
               DO 350 I=1,NWDS                                                                                                      
                  PRES = SIGH(I)*0.01E0*PRESS(I)                                                                                    
                  T    = G(I)                                                                                                       
                  ES   = MAX(F(I),SMIN)                                                                                             
                  IF(ES.GE.SREF) THEN                                                                                               
                    SHUM = ES                                                                                                       
                  ELSE                                                                                                              
                    SHUM = SREF*(EXP(PINV*(1.E0-(SREF/ES)**SPOW)))                                                                  
                  ENDIF                                                                                                             
                  F(I )= SHUM                                                                                                       
                  E    = EXP(A-B/T)                                                                                                 
                  SAT  = EPS1*E/(PRES-EPS2*E)                                                                                       
                  G(I) = SHUM/SAT                                                                                                   
  350          CONTINUE                                                                                                             
                                                                                                                                    
            ELSE                                                                                                                    
               CALL                                XIT('GSHUMH',-9)                                                                 
            ENDIF                                                                                                                   
         ENDIF                                                                                                                      
                                                                                                                                    
C        * PUT SHUM ON FILE 8 AND RHUM ON FILE 9.                                                                                   
                                                                                                                                    
         IBUF(4)=LETAT(L)                                                                                                           
         IBUF(3)=NC4TO8("SHUM")                                                                                                     
         CALL PUTFLD2 (8,F,IBUF,MAXX)                                                                                               
         IF(NRECS.EQ.0) WRITE(6,6025) IBUF                                                                                          
                                                                                                                                    
         IBUF(3)=NC4TO8("RHUM")                                                                                                     
         CALL PUTFLD2 (9,G,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('GSHUMH',-10)                                                                
C---------------------------------------------------------------------                                                              
 5000 FORMAT(10X,2(1X,A4),3E10.0)                                               J4                                                  
 6000 FORMAT(' COORD = ',A4,', MOISTURE VARIABLE = ',A4,                                                                            
     1       ', P.LID (PA)=',E10.3)                                                                                                 
 6005 FORMAT (' SREF=',E12.4,', SPOW=',F10.3)                                                                                       
 6010 FORMAT('0 GSHUMH TRANSFORMED',I6, '  PAIRS')                                                                                  
 6015 FORMAT(' NAME=',A4,',ILEV=',I5,/,' ETA LEVELS=',10I6/(12X,10I6))                                                              
 6025 FORMAT(' ',A4,I10,2X,A4,5I10)                                                                                                 
      END                                                                                                                           
