!> @file
!
!> SUBPROGRAM:    MDL2P       VERT INTRP OF MODEL LVLS TO PRESSURE
!!   PRGRMMR: BLACK           ORG: W/NP22     DATE: 99-09-23       
!!     
!! ABSTRACT:
!!     FOR MOST APPLICATIONS THIS ROUTINE IS THE WORKHORSE
!!     OF THE POST PROCESSOR.  IN A NUTSHELL IT INTERPOLATES
!!     DATA FROM MODEL TO PRESSURE SURFACES.  IT ORIGINATED
!!     FROM THE VERTICAL INTERPOLATION CODE IN THE OLD ETA
!!     POST PROCESSOR SUBROUTINE OUTMAP AND IS A REVISION
!!     OF SUBROUTINE ETA2P.
!!     
!! PROGRAM HISTORY LOG:
!!   99-09-23  T BLACK       - REWRITTEN FROM ETA2P
!!   01-10-25  H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
!!   02-06-12  MIKE BALDWIN - WRF VERSION
!!   02-07-29  H CHUANG - ADD UNDERGROUND FIELDS AND MEMBRANE SLP FOR WRF
!!   04-11-24  H CHUANG - ADD FERRIER'S HYDROMETEOR FIELD
!!   11-02064  J WANG - ADD GRIB2 option
!!   20-03-25  J MENG - remove grib1
!!   21-03-11  B Cui - change local arrays to dimension (im,jsta:jend)
!!   21-10-14  J MENG - 2D DECOMPOSITION
!!   2022-09-01 S Trahan - fixed bugs where extreme atmospheric conditions can cause out-of-bounds access
!!  
!! USAGE:    CALL MDL2P
!!   INPUT ARGUMENT LIST:
!!
!!   OUTPUT ARGUMENT LIST: 
!!     NONE       
!!     
!!   OUTPUT FILES:
!!     NONE
!!     
!!   SUBPROGRAMS CALLED:
!!     UTILITIES:
!!       SCLFLD   - SCALE ARRAY ELEMENTS BY CONSTANT.
!!       CALPOT   - COMPUTE POTENTIAL TEMPERATURE.
!!       CALRH    - COMPUTE RELATIVE HUMIDITY.
!!       CALDWP   - COMPUTE DEWPOINT TEMPERATURE.
!!       BOUND    - BOUND ARRAY ELEMENTS BETWEEN LOWER AND UPPER LIMITS.
!!       CALMCVG  - COMPUTE MOISTURE CONVERGENCE.
!!       CALVOR   - COMPUTE ABSOLUTE VORTICITY.
!!       CALSTRM  - COMPUTE GEOSTROPHIC STREAMFUNCTION.
!!
!!     LIBRARY:
!!       COMMON   - CTLBLK
!!                  RQSTFLD
!!     
!!   ATTRIBUTES:
!!     LANGUAGE: FORTRAN 90
!!     MACHINE : IBM SP
!!
      SUBROUTINE MDL2SIGMA

!
!
      use vrbls3d,       only: pint, t, q, zint, alpint, pmid, exch_h, uh, &
                               vh, omga, q2, cwm, qqw, qqi, qqr, qqs, cfr, &
                               f_rimef, pmidv
!     use vrbls2d,       only:
      use masks,         only: lmh
      use params_mod,    only: d50 , pq0, a2, a3, a4, h1, d01, d608, rgamog,&
                               h1m12, d00, h2, rd, g, gi, h99999
      use ctlblk_mod,    only: jsta_2l, jend_2u, spval, lp1, jsta, jend, lm, &
                               grib, cfld, datapd, fld_info, me, jend_m, im, &
                               jm, im_jm, ista, iend, ista_2l, iend_2u, ista_m, iend_m
      use rqstfld_mod,   only: iget, lvls, id, iavblfld, lvlsxml
      use gridspec_mod,  only :gridtype
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       implicit none
!     
!     INCLUDE MODEL DIMENSIONS.  SET/DERIVE OTHER PARAMETERS.
!     GAMMA AND RGAMOG ARE USED IN THE EXTRAPOLATION OF VIRTUAL
!     TEMPERATURES BEYOND THE UPPER OF LOWER LIMITS OF DATA.
!     
      integer,PARAMETER :: LSIG=22
      real,PARAMETER :: PTSIGO=1.0E4
!     
!     DECLARE VARIABLES.
!     
      LOGICAL READTHK
      LOGICAL IOOMG,IOALL
      LOGICAL DONEFSL1,TSLDONE
      real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: FSL, TSL, QSL, OSL, USL, VSL, Q2SL, &
                                             FSL1, CFRSIG, EGRID1, EGRID2
      REAL GRID1(IM,JM)
      real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: grid2

      REAL SIGO(LSIG+1),DSIGO(LSIG),ASIGO(LSIG)
!
      INTEGER,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: NL1X,NL1XF
!
!
!--- Definition of the following 2D (horizontal) dummy variables
!
!  C1D   - total condensate
!  QW1   - cloud water mixing ratio
!  QI1   - cloud ice mixing ratio
!  QR1   - rain mixing ratio
!  QS1   - snow mixing ratio
!
      real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: C1D, QW1, QI1, QR1, QS1, QG1, AKH
!
      integer I,J,L,LL,LP,LLMH,II,JJ,JJB,JJE,NHOLD
      real PFSIGO,APFSIGO,PSIGO,APSIGO,PNL1,PU,ZU,TU,QU,QSAT,  &
           RHU,TVRU,TVRABV,TABV,QABV,B,AHF,FAC,PL,ZL,TL,QL,    &
           RHL,TMT0,AI,BI,TVRL,TVRBLO,TBLO,QBLO,FACT,          &
           PX,BF,FACF,AHFF,DPSIG,TV,PDV,DENOM,DENOMF,PNL1F,DUM
!
!     
!******************************************************************************
!
!     START MDL2P. 
!     
!     SET TOTAL NUMBER OF POINTS ON OUTPUT GRID.
!
!---------------------------------------------------------------
!
!     *** PART I ***
!
!     VERTICAL INTERPOLATION OF EVERYTHING ELSE.  EXECUTE ONLY
!     IF THERE'S SOMETHING WE WANT.
!
      IF((IGET(205)>0).OR.(IGET(206)>0).OR.              &
         (IGET(207)>0).OR.(IGET(208)>0).OR.              &
         (IGET(209)>0).OR.(IGET(210)>0).OR.              &
         (IGET(216)>0).OR.(IGET(217)>0).OR.              &
         (IGET(211)>0).OR.(IGET(212)>0).OR.              &
         (IGET(213)>0).OR.(IGET(214)>0).OR.              &
         (IGET(215)>0).OR.(IGET(222)>0).OR.              &
         (IGET(243)>0)  ) THEN  !!Air Quality (Plee Oct2003)
!
!---------------------------------------------------------------------
!
!---  VERTICAL INTERPOLATION OF GEOPOTENTIAL, SPECIFIC HUMIDITY, TEMPERATURE, 
!     OMEGA, TKE, & CLOUD FIELDS.  START AT THE UPPERMOST TARGET SIGMA LEVEL.
!
        READTHK=.FALSE.
        IF(READTHK)THEN   ! EITHER READ DSG THICKNESS
	 READ(41)DSIGO  !DSIGO FROM TOP TO BOTTOM
!
         SIGO(1)=0.0
	 DO L=2,LSIG+1
          SIGO(L)=SIGO(L-1)+DSIGO(LSIG-L+2)
	 END DO 
         SIGO(LSIG+1)=1.0
         DO L=1,LSIG
          ASIGO(L)=0.5*(SIGO(L)+SIGO(L+1))
         END DO
        ELSE  ! SPECIFY SIGO
         ASIGO( 1)=   0.0530
         ASIGO( 2)=   0.1580
         ASIGO( 3)=   0.2605
         ASIGO( 4)=   0.3595
         ASIGO( 5)=   0.4550
         ASIGO( 6)=   0.5470
         ASIGO( 7)=   0.6180
         ASIGO( 8)=   0.6690
         ASIGO( 9)=   0.7185
         ASIGO(10)=   0.7585
         ASIGO(11)=   0.7890
         ASIGO(12)=   0.8190
         ASIGO(13)=   0.8480
         ASIGO(14)=   0.8755
         ASIGO(15)=   0.9015
         ASIGO(16)=   0.9260
         ASIGO(17)=   0.9490
         ASIGO(18)=   0.9650
         ASIGO(19)=   0.9745
         ASIGO(20)=   0.9835
         ASIGO(21)=   0.9915
         ASIGO(22)=   0.9975 
!
         SIGO( 1)=   0.0
         SIGO( 2)=   0.1060
         SIGO( 3)=   0.2100
         SIGO( 4)=   0.3110
         SIGO( 5)=   0.4080
         SIGO( 6)=   0.5020
         SIGO( 7)=   0.5920
         SIGO( 8)=   0.6440
         SIGO( 9)=   0.6940
         SIGO(10)=   0.7430
         SIGO(11)=   0.7740
         SIGO(12)=   0.8040
         SIGO(13)=   0.8340
         SIGO(14)=   0.8620
         SIGO(15)=   0.8890
         SIGO(16)=   0.9140
         SIGO(17)=   0.9380
         SIGO(18)=   0.9600
         SIGO(19)=   0.9700
         SIGO(20)=   0.9790
         SIGO(21)=   0.9880
         SIGO(22)=   0.9950 
         SIGO(23)=   1.0
        END IF
! OBTAIN GEOPOTENTIAL AT 1ST LEVEL
       DO J=JSTA_2L,JEND_2U
       DO I=ISTA_2L,IEND_2U
	FSL(I,J)=SPVAL
	AKH(I,J)=SPVAL
        NL1XF(I,J)=LP1
        DO L=1,LP1
        IF(NL1XF(I,J)==LP1.AND.PINT(I,J,L)>PTSIGO)THEN
          NL1XF(I,J)=L
        ENDIF
        ENDDO
       END DO	
       END DO
       DO 167 J=JSTA,JEND
        DO 167 I=ISTA_2L,IEND_2U
	 DONEFSL1=.FALSE.
         PFSIGO=PTSIGO
         APFSIGO=LOG(PFSIGO)
         PNL1=PINT(I,J,NL1XF(I,J))
	 LL=NL1XF(I,J)
         LLMH = NINT(LMH(I,J))
         IF(NL1XF(I,J)==1 .AND. T(I,J,1)<SPVAL           &
      	 .AND. T(I,J,2)<SPVAL .AND. Q(I,J,1)<SPVAL      &
      	 .AND. Q(I,J,2)<SPVAL)THEN
           PU=PINT(I,J,2)
           ZU=ZINT(I,J,2)
           TU=D50*(T(I,J,1)+T(I,J,2))
           QU=D50*(Q(I,J,1)+Q(I,J,2))
           QSAT=PQ0/PU*EXP(A2*(TU-A3)/(TU-A4))
           RHU =QU/QSAT
           IF(RHU>H1)THEN
             RHU=H1
             QU =RHU*QSAT
           ENDIF
           IF(RHU<D01)THEN
             RHU=D01
             QU =RHU*QSAT
           ENDIF
!
           TVRU=TU*(H1+D608*QU)
           TVRABV=TVRU*(PFSIGO/PU)**RGAMOG
           TABV=TVRABV/(H1+D608*QU)
           QSAT=PQ0/PFSIGO*EXP(A2*(TABV-A3)/(TABV-A4))
           QABV =RHU*QSAT
           QABV =MAX(H1M12,QABV)
!==        B    =TABV
           B    =TVRABV       !Marina Tsidulko Dec22, 2003
           AHF  =D00
           FAC  =D00
	   DONEFSL1=.TRUE.
         ELSEIF(NL1XF(I,J)==LP1 .AND. T(I,J,LM-1)<SPVAL     &
      	 .AND. T(I,J,LM-2)<SPVAL .AND. Q(I,J,LM-1)<SPVAL   &
      	 .AND. Q(I,J,LM-2)<SPVAL)THEN
           PL=PINT(I,J,LM-1)
           ZL=ZINT(I,J,LM-1)
           TL=D50*(T(I,J,LM-2)+T(I,J,LM-1))
           QL=D50*(Q(I,J,LM-2)+Q(I,J,LM-1))
!
!--- RH w/r/t water for all conditions  (Ferrier, 25 Jan 02)
!
           QSAT=PQ0/PL*EXP(A2*(TL-A3)/(TL-A4))
           RHL=QL/QSAT
           IF(RHL>H1)THEN
            RHL=H1
            QL =RHL*QSAT
           ENDIF
           IF(RHL<D01)THEN
             RHL=D01
             QL =RHL*QSAT
           ENDIF
!
           TVRL  =TL*(H1+D608*QL)
           TVRBLO=TVRL*(PFSIGO/PL)**RGAMOG
           TBLO  =TVRBLO/(H1+D608*QL)
           QSAT=PQ0/PFSIGO*EXP(A2*(TBLO-A3)/(TBLO-A4))
           QBLO =RHL*QSAT
           QBLO =MAX(H1M12,QBLO)
!==        B    =TBLO
           B    =TVRBLO             !Marina Tsidulko Dec22, 2003
           AHF  =D00
           FAC  =D00
	   DONEFSL1=.TRUE.
         ELSEIF(NL1XF(I,J)<LP1) THEN ! must check the "<LP1" first; do not combine these into a single "if"
          IF(T(I,J,NL1XF(I,J))<SPVAL                 &
     &	   .AND. Q(I,J,NL1XF(I,J))<SPVAL)THEN
!==        B     =T(I,J,NL1XF(I,J)) !Marina Tsidulko Dec22, 2003
           B     =T(I,J,NL1XF(I,J))*(H1+D608*Q(I,J,NL1XF(I,J)))
           DENOM=(ALPINT(I,J,NL1XF(I,J)+1)-ALPINT(I,J,NL1XF(I,J)-1))
!==        AHF  =(B-T(I,J,NL1XF(I,J)-1))/DENOM
           AHF  =(B-T(I,J,NL1XF(I,J)-1)*(H1+D608*Q(I,J,NL1XF(I,J)-1))) &
                   /DENOM            !Marina Tsidulko Dec22, 2003
           FAC =H2*LOG(PMID(I,J,NL1XF(I,J)))
	   DONEFSL1=.TRUE.
          END IF
         END IF
!
         if(DONEFSL1)FSL1(I,J)=(PNL1-PFSIGO)/(PFSIGO+PNL1)              &
                 *((APFSIGO+ALPINT(I,J,NL1XF(I,J))-FAC)*AHF+B)*RD*H2    &
                 +ZINT(I,J,NL1XF(I,J))*G
! COMPUTE EXCHANGE COEFFICIENT ON FIRST INTERFACTE
         IF(NL1XF(I,J)<=2 .OR. NL1XF(I,J)>(LLMH+1))THEN 
          AKH(I,J)=0.0
         ELSE
          FACT=(APFSIGO-LOG(PINT(I,J,LL)))/                            &
     &         (LOG(PINT(I,J,LL))-LOG(PINT(I,J,LL-1)))
! EXCH_H is on the bottom of model interfaces
          IF(EXCH_H(I,J,LL-2)<SPVAL .AND. EXCH_H(I,J,LL-1)<SPVAL) &
     &      AKH(I,J)=EXCH_H(I,J,LL-1)+(EXCH_H(I,J,LL-1)                 &
     &	    -EXCH_H(I,J,LL-2))*FACT 
	 END IF    
 167    CONTINUE
! OUTPUT FIRST LAYER GEOPOTENTIAL
!       GEOPOTENTIAL (SCALE BY GI)
        IF (IGET(205)>0) THEN
         IF (LVLS(1,IGET(205))>0) THEN
!$omp  parallel do
           DO J=JSTA,JEND
           DO I=ISTA,IEND
             IF(FSL1(I,J)<SPVAL) THEN
                GRID1(I,J)=FSL1(I,J)*GI
             ELSE
                GRID1(I,J)=SPVAL
             ENDIF
           ENDDO
           ENDDO
          if(grib=='grib2')then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(205))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
         ENDIF
        ENDIF

! OUTPUT FIRST INTERFACE KH Heat Diffusivity
         IF (IGET(243)>0) THEN  !!Air Quality (Plee Oct2003) ^^^^^
          IF (LVLS(1,IGET(243))>0) THEN
!$omp  parallel do
           DO J=JSTA,JEND
           DO I=ISTA,IEND
             GRID1(I,J)=AKH(I,J)
           ENDDO
           ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(243))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif

           if(me==0)print*,'output Heat Diffusivity'
         ENDIF
        ENDIF 	

!***
!***  BECAUSE SIGMA LAYERS DO NOT GO UNDERGROUND,  DO ALL
!***  INTERPOLATION ABOVE GROUND NOW.
!***
!
        DO 310 LP=1,LSIG
        NHOLD=0
!
        DO J=JSTA_2L,JEND_2U
        DO I=ISTA_2L,IEND_2U

!
        TSL(I,J)=SPVAL
        QSL(I,J)=SPVAL
        FSL(I,J)=SPVAL
        OSL(I,J)=SPVAL
        USL(I,J)=SPVAL
        VSL(I,J)=SPVAL
        Q2SL(I,J)=SPVAL
        C1D(I,J)=SPVAL      ! Total condensate
        QW1(I,J)=SPVAL      ! Cloud water
        QI1(I,J)=SPVAL      ! Cloud ice
        QR1(I,J)=SPVAL      ! Rain 
        QS1(I,J)=SPVAL      ! Snow (precip ice) 
	QG1(I,J)=SPVAL
        CFRSIG(I,J)=SPVAL
!
!***  LOCATE VERTICAL INDEX OF MODEL MIDLAYER JUST BELOW
!***  THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
!
        NL1X(I,J)=LP1
        DO L=2,LM
        LLMH = NINT(LMH(I,J))
        PSIGO=PTSIGO+ASIGO(LP)*(PINT(I,J,LLMH+1)-PTSIGO)
        IF(NL1X(I,J)==LP1.AND.PMID(I,J,L)>PSIGO)THEN
          NL1X(I,J)=L
        ENDIF
        ENDDO
!
!  IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
!  BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
!  WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
!  WILL EXTRAPOLATE TO THAT POINT
!
        IF(NL1X(I,J)==LP1.AND.PINT(I,J,LLMH+1)>=PSIGO)THEN
          NL1X(I,J)=LM
        ENDIF
!
!        if(NL1X(I,J)==LP1)print*,'Debug: NL1X=LP1 AT '
!     1 ,i,j,lp
        ENDDO
        ENDDO
!
!mptest        IF(NHOLD==0)GO TO 310
!
!$omp  parallel do private(i,j,ll,llmh,psigo,apsigo,fact,dum,pl,     &
!$omp &         zl,tl,ql,ai,bi,qsat,rhl,tvrl,tvrblo,tblo,tmt0,       &
!$omp &         qblo,pnl1,fac,ahf)
!hc        DO 220 NN=1,NHOLD
!hc        I=IHOLD(NN)
!hc        J=JHOLD(NN)
        DO 220 J=JSTA,JEND      ! Moorthi on Nov 26 2014
!       DO 220 J=JSTA_2L,JEND_2U
        DO 220 I=ISTA,IEND
        LL=NL1X(I,J)
!---------------------------------------------------------------------
!***  VERTICAL INTERPOLATION OF GEOPOTENTIAL, TEMPERATURE, SPECIFIC
!***  HUMIDITY, CLOUD WATER/ICE, OMEGA, WINDS, AND TKE.
!---------------------------------------------------------------------
!
!HC        IF(NL1X(I,J)<=LM)THEN
        LLMH = NINT(LMH(I,J))
	PSIGO=PTSIGO+ASIGO(LP)*(PINT(I,J,LLMH+1)-PTSIGO) 
	APSIGO=LOG(PSIGO)
        IF(NL1X(I,J)<=LLMH)THEN
!
!---------------------------------------------------------------------
!          INTERPOLATE LINEARLY IN LOG(P)
!***  EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
!***  INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
!***  EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
!---------------------------------------------------------------------
!

          FACT=(APSIGO-LOG(PMID(I,J,LL)))/                         &
     &         (LOG(PMID(I,J,LL))-LOG(PMID(I,J,LL-1)))
          TSL(I,J)=T(I,J,LL)+(T(I,J,LL)-T(I,J,LL-1))*FACT
          IF(Q(I,J,LL)<SPVAL .AND. Q(I,J,LL-1)<SPVAL)         &
     &      QSL(I,J)=Q(I,J,LL)+(Q(I,J,LL)-Q(I,J,LL-1))*FACT
	  IF(gridtype=='A')THEN
           IF(UH(I,J,LL)<SPVAL .AND. UH(I,J,LL-1)<SPVAL)      &
     &       USL(I,J)=UH(I,J,LL)+(UH(I,J,LL)-UH(I,J,LL-1))*FACT           
           IF(VH(I,J,LL)<SPVAL .AND. VH(I,J,LL-1)<SPVAL)      &
     &       VSL(I,J)=VH(I,J,LL)+(VH(I,J,LL)-VH(I,J,LL-1))*FACT
	  END IF 
          IF(OMGA(I,J,LL)<SPVAL .AND. OMGA(I,J,LL-1)<SPVAL)   &
     &      OSL(I,J)=OMGA(I,J,LL)+(OMGA(I,J,LL)-OMGA(I,J,LL-1))*FACT
          IF(Q2(I,J,LL)<SPVAL .AND. Q2(I,J,LL-1)<SPVAL)       &
     &      Q2SL(I,J)=Q2(I,J,LL)+(Q2(I,J,LL)-Q2(I,J,LL-1))*FACT      
!          FSL(I,J)=ZMID(I,J,LL)+(ZMID(I,J,LL)-ZMID(I,J,LL-1))*FACT
!          FSL(I,J)=FSL(I,J)*G
!hc          QSAT=PQ0/PSIGO
!hc     1      *EXP(A2*(TSL(I,J)-A3)/(TSL(I,J)-A4))
!
!hc          RHL=QSL(I,J)/QSAT
!
!hc          IF(RHL>1.) QSL(I,J)=QSAT
!hc          IF(RHL<0.01) QSL(I,J)=0.01*QSAT
          IF(Q2SL(I,J)<0.0) Q2SL(I,J)=0.0
!	  
!HC ADD FERRIER'S HYDROMETEOR
          IF(CWM(I,J,LL)<SPVAL .AND. CWM(I,J,LL-1)<SPVAL)    &
     &       C1D(I,J)=CWM(I,J,LL)+(CWM(I,J,LL)-CWM(I,J,LL-1))*FACT          
          C1D(I,J)=MAX(C1D(I,J),H1M12)      ! Total condensate	  
          IF(QQW(I,J,LL)<SPVAL .AND. QQW(I,J,LL-1)<SPVAL)    &
     &       QW1(I,J)=QQW(I,J,LL)+(QQW(I,J,LL)-QQW(I,J,LL-1))*FACT
          QW1(I,J)=MAX(QW1(I,J),H1M12)      ! Cloud water
          IF(QQI(I,J,LL)<SPVAL .AND. QQI(I,J,LL-1)<SPVAL)    &
     &       QI1(I,J)=QQI(I,J,LL)+(QQI(I,J,LL)-QQI(I,J,LL-1))*FACT
          QI1(I,J)=MAX(QI1(I,J),H1M12)      ! Cloud ice
          IF(QQR(I,J,LL)<SPVAL .AND. QQR(I,J,LL-1)<SPVAL)    &
     &       QR1(I,J)=QQR(I,J,LL)+(QQR(I,J,LL)-QQR(I,J,LL-1))*FACT
          QR1(I,J)=MAX(QR1(I,J),H1M12)      ! Rain 
          IF(QQS(I,J,LL)<SPVAL .AND. QQS(I,J,LL-1)<SPVAL)    &
     &       QS1(I,J)=QQS(I,J,LL)+(QQS(I,J,LL)-QQS(I,J,LL-1))*FACT
          QS1(I,J)=MAX(QS1(I,J),H1M12)      ! Snow (precip ice) 
          IF(CFR(I,J,LL)<SPVAL .AND. CFR(I,J,LL-1)<SPVAL)    &
     &       CFRSIG(I,J)=CFR(I,J,LL)+(CFR(I,J,LL)-CFR(I,J,LL-1))*FACT
          CFRSIG(I,J)=MAX(CFRSIG(I,J),H1M12)
          IF(QQS(I,J,LL)<SPVAL .AND. QQS(I,J,LL-1)<SPVAL)THEN
	   DUM=F_RimeF(I,J,LL)+(F_RimeF(I,J,LL)-F_RimeF(I,J,LL-1))*FACT
           IF(DUM <= 5.0)THEN
	    QG1(I,J)=H1M12
	   ELSE
	    QG1(I,J)=QS1(I,J)
	    QS1(I,J)=H1M12
	   END IF 
	  END IF  
! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE 
! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
! GOUND
        ELSE
!          ii=91
!          jj=13
!          if(i==ii.and.j==jj)print*,'Debug: underg extra at i,j,lp' &
!     &,   i,j,lp
	  PL=PINT(I,J,LM-1)
          ZL=ZINT(I,J,LM-1)
          TL=0.5*(T(I,J,LM-2)+T(I,J,LM-1))
          QL=0.5*(Q(I,J,LM-2)+Q(I,J,LM-1))
          AI=0.008855
          BI=1.
          TMT0=TL-A3
          IF(TMT0<-20.)THEN
            AI=0.007225
            BI=0.9674
          ENDIF
          QSAT=PQ0/PL*EXP(A2*(TL-A3)/(TL-A4))
!
          RHL=QL/QSAT
!
          IF(RHL>1.)THEN
            RHL=1.
            QL =RHL*QSAT
          ENDIF
!
          IF(RHL<0.01)THEN
            RHL=0.01
            QL =RHL*QSAT
          ENDIF
!
          TVRL  =TL*(1.+0.608*QL)
          TVRBLO=TVRL*(PSIGO/PL)**RGAMOG
          TBLO  =TVRBLO/(1.+0.608*QL)
!     
          TMT0=TBLO-A3
          AI=0.008855
          BI=1.
          IF(TMT0<-20.)THEN
            AI=0.007225
            BI=0.9674
          ENDIF
          QSAT=PQ0/PSIGO*EXP(A2*(TBLO-A3)/(TBLO-A4))
!
!          TSL(I,J)=TBLO
          QBLO     = RHL*QSAT
          QSL(I,J) = MAX(1.E-12,QBLO)
          IF(gridtype=='A')THEN
            USL(I,J) = UH(I,J,LLMH)
            VSL(I,J) = VH(I,J,LLMH)
          END IF 
          OSL(I,J)  = OMGA(I,J,LLMH)
          Q2SL(I,J) = max(0.0,0.5*(Q2(I,J,LLMH-1)+Q2(I,J,LLMH)))
          PNL1      = PINT(I,J,NL1X(I,J))
          FAC       = 0.
          AHF       = 0.0
!
!--- Set hydrometeor fields to zero below ground
          C1D(I,J)=0.
          QW1(I,J)=0.
          QI1(I,J)=0.
          QR1(I,J)=0.
          QS1(I,J)=0.
	  QG1(I,J)=0.
          CFRSIG(I,J)=0.
        END IF
  220   CONTINUE
!
! OBTAIN GEOPOTENTIAL AND KH ON INTERFACES 
       DO J=JSTA_2L,JEND_2U
       DO I=ISTA_2L,IEND_2U
        FSL(I,J)=SPVAL
	AKH(I,J)=SPVAL
        NL1XF(I,J)=LP1
        LLMH = NINT(LMH(I,J))
        PSIGO=PTSIGO+SIGO(LP+1)*(PINT(I,J,LLMH+1)-PTSIGO)
        DO L=1,LP1
        IF(NL1XF(I,J)==LP1.AND.PINT(I,J,L)>PSIGO)THEN
          NL1XF(I,J)=L
        ENDIF
        ENDDO
       END DO
       END DO
!
!      DO J=JSTA_2L,JEND_2U
       DO J=JSTA,JEND          ! Moorthi on 26 Nov 2014
       DO I=ISTA,IEND
          DONEFSL1=.FALSE.
	  TSLDONE=.FALSE.
          LLMH = NINT(LMH(I,J))  
          PFSIGO=PTSIGO+SIGO(LP+1)*(PINT(I,J,LLMH+1)-PTSIGO) 
          PSIGO=PTSIGO+ASIGO(LP)*(PINT(I,J,LLMH+1)-PTSIGO)
          APFSIGO=LOG(PFSIGO)
          PNL1F=PINT(I,J,NL1XF(I,J))   
	  LL=NL1XF(I,J)
          IF(NL1XF(I,J)==1 .AND. T(I,J,1)<SPVAL        &
     &	 .AND. T(I,J,2)<SPVAL .AND. Q(I,J,1)<SPVAL    &
     &	 .AND. Q(I,J,2)<SPVAL)THEN
              PU=PINT(I,J,2)
              ZU=ZINT(I,J,2)
              TU=D50*(T(I,J,1)+T(I,J,2))
              QU=D50*(Q(I,J,1)+Q(I,J,2))
!
!--- RH w/r/t water for all conditions
!
              QSAT=PQ0/PU*EXP(A2*(TU-A3)/(TU-A4))
              RHU =QU/QSAT
              IF(RHU>H1)THEN
                RHU=H1
                QU =RHU*QSAT
              ENDIF
              IF(RHU<D01)THEN
                RHU=D01
                QU =RHU*QSAT
              ENDIF
!
              TVRU=TU*(H1+D608*QU)
!              TVRABV=TVRU*(SPL(L)/PU)**RGAMOG
!==           TVRABV=TVRU*(PFSIGO/PU)**RGAMOG
              PX=(PFSIGO+PNL1F)*0.5      !Marina Tsidulko Dec22, 2003
              TVRABV=TVRU*(PX/PU)**RGAMOG!Marina Tsidulko Dec22, 2003
              TABV=TVRABV/(H1+D608*QU)
! ADD ADDITIONAL BLOCK FOR INTERPOLATING HEIGHT ONTO INTERFACES
!==           BF   =TABV
              BF   =TVRABV     !Marina Tsidulko Dec22, 2003
              FACF =D00
              AHFF =D00
	      DONEFSL1=.TRUE.
!
!
            ELSEIF(NL1XF(I,J)==LP1 .AND. T(I,J,LM-1)<SPVAL    &
     &	 .AND. T(I,J,LM-2)<SPVAL .AND. Q(I,J,LM-1)<SPVAL     &
     &	 .AND. Q(I,J,LM-2)<SPVAL)THEN
!
!           EXTRAPOLATION AT LOWER BOUND.  THE LOWER BOUND IS
!           LM IF OLDRD=.FALSE.  IF OLDRD=.TRUE. THE LOWER
!           BOUND IS THE FIRST ATMOSPHERIC ETA LAYER.
!
              PL=PINT(I,J,LM-1)
              ZL=ZINT(I,J,LM-1)
              TL=D50*(T(I,J,LM-2)+T(I,J,LM-1))
              QL=D50*(Q(I,J,LM-2)+Q(I,J,LM-1)) 
!
!--- RH w/r/t water for all conditions  (Ferrier, 25 Jan 02)
!
              QSAT=PQ0/PL*EXP(A2*(TL-A3)/(TL-A4))
              RHL=QL/QSAT
              IF(RHL>H1)THEN
               RHL=H1
               QL =RHL*QSAT
              ENDIF
              IF(RHL<D01)THEN
                RHL=D01
                QL =RHL*QSAT
              ENDIF
!
              TVRL  =TL*(H1+D608*QL)
!==           TVRBLO=TVRL*(PFSIGO/PL)**RGAMOG
              PX=(PFSIGO+PNL1F)*0.5       !Marina Tsidulko Dec22, 2003
              TVRBLO=TVRL*(PX/PL)**RGAMOG !Marina Tsidulko Dec22, 2003
              TBLO  =TVRBLO/(H1+D608*QL)
! ADD ADDITIONAL BLOCK FOR INTERPOLATING HEIGHT ONTO INTERFACES
!==           BF   =TBLO
              BF   =TVRBLO    !Marina Tsidulko Dec22, 2003
              FACF =D00
              AHFF =D00
	      DONEFSL1=.TRUE.
	      TSLDONE=.TRUE.
!
!
            ELSEIF(NL1XF(I,J)<LP1) THEN ! must check the "<LP1" first; do not combine these into a single "if"
             IF(T(I,J,NL1XF(I,J))<SPVAL              &
     &	  .AND. Q(I,J,NL1XF(I,J))<SPVAL)THEN
!
!           INTERPOLATION BETWEEN LOWER AND UPPER BOUNDS.
!
! ADD ADDITIONAL BLOCK FOR INTERPOLATING HEIGHT ONTO INTERFACES
!             if(NL1XF(I,J)==LMp1)
!    +        print*,'Debug: using bad temp at',i,j
!==           BF   =T(I,J,NL1XF(I,J))  !Marina Tsidulko Dec22, 2003
              BF   =T(I,J,NL1XF(I,J))*(H1+D608*Q(I,J,NL1XF(I,J)))
!HC              FACF =H2*LOG(PT+PDSL(I,J)*AETA(NL1XF(I,J)))
              FACF =H2*LOG(PMID(I,J,NL1XF(I,J)))
              DENOMF=(ALPINT(I,J,NL1XF(I,J)+1)-ALPINT(I,J,NL1XF(I,J)-1))
!==           AHFF=(BF-T(I,J,NL1XF(I,J)-1))/DENOMF
           AHFF=(BF-T(I,J,NL1XF(I,J)-1)*(H1+D608*Q(I,J,NL1XF(I,J)-1)))  &
                   /DENOMF             !Marina Tsidulko Dec22, 2003
!
              DONEFSL1=.TRUE.     
             ENDIF
            ENDIF   
!
            IF(DONEFSL1)then
	     FSL(I,J)=(PNL1F-PFSIGO)/(PFSIGO+PNL1F)                   &
              *((APFSIGO+ALPINT(I,J,NL1XF(I,J))-FACF)*AHFF+BF)*RD*H2  &
                 +ZINT(I,J,NL1XF(I,J))*G
!  OBTAIN TEMPERATURE AT MID-SIGMA LAYER BASED ON HYDROSTATIC
             DPSIG=(SIGO(LP+1)-SIGO(LP))*(PINT(I,J,LLMH+1)-PTSIGO)
!           IF(J==jj.and.i==ii)print*,'Debug:L,OLD T= ',
!    +      L,TSL(I,J)
             IF(.NOT.TSLDONE) THEN
               TSL(I,J)=(FSL1(I,J)-FSL(I,J))*PSIGO/(RD*DPSIG)
             ENDIF
             FSL1(I,J)=FSL(I,J)
	     IF(.NOT.TSLDONE) THEN
               TV=TSL(I,J)
               TSL(I,J)=TV/(H1+D608*QSL(I,J))
             ENDIF
             QSAT=PQ0/PSIGO *EXP(A2*(TSL(I,J)-A3)/(TSL(I,J)-A4))
!
             RHL=QSL(I,J)/QSAT
!
             IF(RHL>1.) QSL(I,J)=QSAT
             IF(RHL<0.01) QSL(I,J)=0.01*QSAT
            END IF
!           IF(J==jj.and.i==ii)print*,'Debug:L,T,Q,Q2,FSL=',
!    +      L,TSL(I,J),QSL(I,J),Q2SL(I,J),FSL(I,J)

! COMPUTE EXCHANGE COEFFICIENT ON INTERFACES
           IF(NL1XF(I,J)<=2 .OR. NL1XF(I,J)>(LLMH+1))THEN
             AKH(I,J)=0.0
           ELSE
            FACT=(APFSIGO-LOG(PINT(I,J,LL)))/             &
     &        (LOG(PINT(I,J,LL))-LOG(PINT(I,J,LL-1)))
! EXCH_H is on the bottom of model interfaces
            IF(EXCH_H(I,J,LL-2)<SPVAL .AND.             &
     &	    EXCH_H(I,J,LL-1)<SPVAL)                     &
     &        AKH(I,J)=EXCH_H(I,J,LL-1)+(EXCH_H(I,J,LL-1)  &
     &         -EXCH_H(I,J,LL-2))*FACT
           END IF

! END OF HYDROSTATIC TEMPERATURE DERIVATION 
       END DO
       END DO
!	
! VERTICAL INTERPOLATION FOR WIND FOR E and B GRIDS
!
        if(gridtype=='B' .or. gridtype=='E') &
	 call exch(PINT(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LP1))  
        IF(gridtype=='E')THEN
        DO J=JSTA,JEND
!        DO I=1,IM-MOD(J,2)
        DO I=ISTA,IEND-MOD(J,2) !Jesse 20211014
!
!***  LOCATE VERTICAL INDEX OF MODEL MIDLAYER FOR V POINT JUST BELOW
!***  THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
!
	 LLMH = NINT(LMH(I,J))

!Jesse 20211014
!         IF(J == 1 .AND. I < IM)THEN   !SOUTHERN BC
!           PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
!         ELSE IF(J==JM .AND. I<IM)THEN   !NORTHERN BC
!           PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
!         ELSE IF(I == 1 .AND. MOD(J,2) == 0) THEN   !WESTERN EVEN BC
!           PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
!         ELSE IF(I == IM .AND. MOD(J,2) == 0) THEN   !EASTERN EVEN BC
!           PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
!         ELSE IF (MOD(J,2) < 1) THEN
!           PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I-1,J,LLMH+1)              &
!     &       +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
!         ELSE
!           PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1)              &
!     &       +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
!         END IF 

         IF(J == 1 .AND. I < IEND)THEN   !SOUTHERN BC
           PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
         ELSE IF(J==JM .AND. I<IEND)THEN   !NORTHERN BC
           PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
         ELSE IF(I == ISTA .AND. MOD(J,2) == 0) THEN   !WESTERN EVEN BC
           PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
         ELSE IF(I == IEND .AND. MOD(J,2) == 0) THEN   !EASTERN EVEN BC
           PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
         ELSE IF (MOD(J,2) < 1) THEN
           PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I-1,J,LLMH+1)              &
     &       +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
         ELSE
           PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1)              &
     &       +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
         END IF

         PSIGO=PTSIGO+ASIGO(LP)*(PDV-PTSIGO)
	 NL1X(I,J)=LP1 
         DO L=2,LM
          IF(NL1X(I,J)==LP1.AND.PMIDV(I,J,L)>PSIGO)THEN
           NL1X(I,J)=L
          ENDIF
         ENDDO
!
!  IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
!  BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
!  WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
!  WILL EXTRAPOLATE TO THAT POINT
!
         IF(NL1X(I,J)==LP1.AND. PDV>PSIGO)THEN	
          NL1X(I,J)=LM
         ENDIF
!
        ENDDO
        ENDDO
!
        DO 230 J=JSTA,JEND
!        DO 230 I=1,IM-MOD(j,2)
        DO 230 I=ISTA,IEND-MOD(j,2) !Jesse 20211014

         LLMH = NINT(LMH(I,J))

!Jesse 20211014
!         IF(J == 1 .AND. I < IM)THEN   !SOUTHERN BC
!           PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
!         ELSE IF(J==JM .AND. I<IM)THEN   !NORTHERN BC
!           PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
!         ELSE IF(I == 1 .AND. MOD(J,2) == 0) THEN   !WESTERN EVEN BC
!           PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
!         ELSE IF(I == IM .AND. MOD(J,2) == 0) THEN   !EASTERN EVEN BC
!           PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
!         ELSE IF (MOD(J,2) < 1) THEN
!           PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I-1,J,LLMH+1)        &
!     &       +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
!         ELSE
!           PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1)        &
!     &       +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
!         END IF

         IF(J == 1 .AND. I < IEND)THEN   !SOUTHERN BC
           PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
         ELSE IF(J==JM .AND. I<IEND)THEN   !NORTHERN BC
           PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
         ELSE IF(I == ISTA .AND. MOD(J,2) == 0) THEN   !WESTERN EVEN BC
           PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
         ELSE IF(I == IEND .AND. MOD(J,2) == 0) THEN   !EASTERN EVEN BC
           PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
         ELSE IF (MOD(J,2) < 1) THEN
           PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I-1,J,LLMH+1)        &
     &       +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
         ELSE
           PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1)        &
     &       +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
         END IF

         PSIGO=PTSIGO+ASIGO(LP)*(PDV-PTSIGO)
	 APSIGO=LOG(PSIGO)  
         LL=NL1X(I,J)
!---------------------------------------------------------------------
!***  VERTICAL INTERPOLATION OF WINDS FOR A-E GRID
!---------------------------------------------------------------------
!         
!HC        IF(NL1X(I,J)<=LM)THEN
         LLMH = NINT(LMH(I,J))
         IF(NL1X(I,J)<=LLMH)THEN
!
!---------------------------------------------------------------------
!          INTERPOLATE LINEARLY IN LOG(P)
!***  EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
!***  INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
!***  EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
!---------------------------------------------------------------------
!
	   
          FACT=(APSIGO-LOG(PMIDV(I,J,LL)))/                   &
     &         (LOG(PMIDV(I,J,LL))-LOG(PMIDV(I,J,LL-1)))
          IF(UH(I,J,LL)<SPVAL .AND. UH(I,J,LL-1)<SPVAL)  &
     &       USL(I,J)=UH(I,J,LL)+(UH(I,J,LL)-UH(I,J,LL-1))*FACT
          IF(VH(I,J,LL)<SPVAL .AND. VH(I,J,LL-1)<SPVAL)  &
     &       VSL(I,J)=VH(I,J,LL)+(VH(I,J,LL)-VH(I,J,LL-1))*FACT
!
! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE 
! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
! GOUND
        ELSE
          IF(UH(I,J,LLMH)<SPVAL)USL(I,J)=UH(I,J,LLMH)
	  IF(VH(I,J,LLMH)<SPVAL)VSL(I,J)=VH(I,J,LLMH)
        END IF
  230   CONTINUE
        JJB=JSTA 
        IF(MOD(JSTA,2)==0)JJB=JSTA+1
        JJE=JEND
        IF(MOD(JEND,2)==0)JJE=JEND-1
        DO J=JJB,JJE,2 !chc
          USL(IEND,J)=USL(IEND-1,J)
	  VSL(IEND,J)=VSL(IEND-1,J)
        END DO
	
	ELSE IF (gridtype=='B')THEN
         DO J=JSTA,JEND_M
         DO I=ISTA,IEND_M
!
!***  LOCATE VERTICAL INDEX OF MODEL MIDLAYER FOR V POINT JUST BELOW
!***  THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
!
           PDV=0.25*(PINT(I,J,LP1)+PINT(I+1,J,LP1)                       &
             +PINT(I,J+1,LP1)+PINT(I+1,J+1,LP1))  

           PSIGO=PTSIGO+ASIGO(LP)*(PDV-PTSIGO)
	   NL1X(I,J)=LP1 
           DO L=2,LM
            IF(NL1X(I,J)==LP1.AND.PMIDV(I,J,L)>PSIGO)THEN
             NL1X(I,J)=L
            ENDIF
           ENDDO
!
!  IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
!  BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
!  WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
!  WILL EXTRAPOLATE TO THAT POINT
!
           IF(NL1X(I,J)==LP1.AND. PDV>PSIGO)THEN	
             NL1X(I,J)=LM
           ENDIF
!
         ENDDO
         ENDDO
!
         DO 231 J=JSTA,JEND_M
         DO 231 I=ISTA,IEND_M
	  PDV=0.25*(PINT(I,J,LP1)+PINT(I+1,J,LP1)                       &
             +PINT(I,J+1,LP1)+PINT(I+1,J+1,LP1))
          PSIGO=PTSIGO+ASIGO(LP)*(PDV-PTSIGO)
	  APSIGO=LOG(PSIGO)  
          LL=NL1X(I,J)
!---------------------------------------------------------------------
!***  VERTICAL INTERPOLATION OF WINDS FOR A-E GRID
!---------------------------------------------------------------------
!         
!HC        IF(NL1X(I,J)<=LM)THEN
          LLMH = NINT(LMH(I,J))
          IF(NL1X(I,J)<=LLMH)THEN
!
!---------------------------------------------------------------------
!          INTERPOLATE LINEARLY IN LOG(P)
!***  EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
!***  INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
!***  EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
!---------------------------------------------------------------------
!
	   
           FACT=(APSIGO-LOG(PMIDV(I,J,LL)))/                   &
     &         (LOG(PMIDV(I,J,LL))-LOG(PMIDV(I,J,LL-1)))
           IF(UH(I,J,LL)<SPVAL .AND. UH(I,J,LL-1)<SPVAL)  &
     &        USL(I,J)=UH(I,J,LL)+(UH(I,J,LL)-UH(I,J,LL-1))*FACT
           IF(VH(I,J,LL)<SPVAL .AND. VH(I,J,LL-1)<SPVAL)  &
     &        VSL(I,J)=VH(I,J,LL)+(VH(I,J,LL)-VH(I,J,LL-1))*FACT
!
! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE 
! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
! GOUND
          ELSE
           IF(UH(I,J,LLMH)<SPVAL)USL(I,J)=UH(I,J,LLMH)
	   IF(VH(I,J,LLMH)<SPVAL)VSL(I,J)=VH(I,J,LLMH)
          END IF
  231   CONTINUE
	
	
	
        END IF  ! END OF WIND INTERPOLATION FOR NMM


!
!---------------------------------------------------------------------
!        LOAD GEOPOTENTIAL AND TEMPERATURE INTO STANDARD LEVEL 
!        ARRAYS FOR THE NEXT PASS.
!---------------------------------------------------------------------
!     
!---------------------------------------------------------------------
!***  CALCULATE 1000MB GEOPOTENTIALS CONSISTENT WITH SLP OBTAINED 
!***  FROM THE MESINGER OR NWS SHUELL SLP REDUCTION.
!---------------------------------------------------------------------
!     
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!        *** PART II ***
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!
!        INTERPOLATE/OUTPUT SELECTED FIELDS.
!
!---------------------------------------------------------------------
!
!***  OUTPUT GEOPOTENTIAL (SCALE BY GI)
!
        IF(IGET(205)>0)THEN
          IF(LVLS(LP+1,IGET(205))>0)THEN
!$omp  parallel do
            DO J=JSTA,JEND
            DO I=ISTA,IEND
              IF(FSL(I,J)<SPVAL) THEN
                GRID1(I,J)=FSL(I,J)*GI
              ELSE
                GRID1(I,J)=SPVAL
              ENDIF
            ENDDO
            ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(205))
            fld_info(cfld)%lvl=LVLSXML(LP+1,IGET(205))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
        ENDIF
!	
!***  OUTPUT EXCHANGE COEEFICIENT
!
! OUTPUT FIRST INTERFACE KH Heat Diffusivity
         IF (IGET(243)>0) THEN  !!Air Quality (Plee Oct2003) ^^^^^
          IF (LVLS(LP+1,IGET(243))>0) THEN
!$omp  parallel do
           DO J=JSTA,JEND
           DO I=ISTA,IEND
             GRID1(I,J)=AKH(I,J)
	     IF(LP==(LSIG+1))GRID1(I,J)=0.0  !! NO SLIP ASSUMTION FOR CMAQ
           ENDDO
           ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(243))
            fld_info(cfld)%lvl=LVLSXML(LP+1,IGET(243))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
           if(me==0)print*,'output Heat Diffusivity'
         ENDIF
        ENDIF 
!     
!***  TEMPERATURE
!
        IF(IGET(206)>0) THEN
          IF(LVLS(LP,IGET(206))>0) THEN
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=TSL(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(206))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(206))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
        ENDIF
!     
!***  PRESSURE
!
        IF(IGET(216)>0)THEN
          IF(LVLS(LP,IGET(216))>0)THEN
!$omp  parallel do
             DO J=JSTA,JEND
             DO I=ISTA,IEND
	       LLMH = NINT(LMH(I,J))  
               GRID1(I,J)=PTSIGO+ASIGO(LP)*(PINT(I,J,LLMH+1)-PTSIGO)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(216))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(216))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
        ENDIF
!
!***  SPECIFIC HUMIDITY.
!
        IF(IGET(207)>0)THEN
          IF(LVLS(LP,IGET(207))>0)THEN
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=QSL(I,J)
             ENDDO
             ENDDO
             CALL BOUND(GRID1,H1M12,H99999)
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(207))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(207))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
        ENDIF
!     
!***  OMEGA
!
        IF(IGET(210)>0)THEN
          IF(LVLS(LP,IGET(210))>0)THEN
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=OSL(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(210))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(210))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
        ENDIF
!     
!***  U AND/OR V WIND
!
        IF(IGET(208)>0.OR.IGET(209)>0)THEN
          IF(LVLS(LP,IGET(208))>0.OR.LVLS(LP,IGET(209))>0) then
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=USL(I,J)
               GRID2(I,J)=VSL(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(208))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(208))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(209))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(209))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID2(ista:iend,jsta:jend)
          endif
          ENDIF
        ENDIF
!     
!***  TURBULENT KINETIC ENERGY
!
         IF (IGET(217)>0) THEN
          IF (LVLS(LP,IGET(217))>0) THEN
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=Q2SL(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(217))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(217))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
         ENDIF
!     
!***  CLOUD WATER
!
         IF (IGET(211)>0) THEN
          IF (LVLS(LP,IGET(211))>0) THEN
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=QW1(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(211))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(211))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
         ENDIF
!
!***  CLOUD ICE 
!
         IF (IGET(212)>0) THEN
          IF (LVLS(LP,IGET(212))>0) THEN
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=QI1(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(212))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(212))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
         ENDIF
!
!---  RAIN
         IF (IGET(213)>0) THEN
          IF (LVLS(LP,IGET(213))>0) THEN 
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=QR1(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(213))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(213))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
         ENDIF
!
!---  SNOW
         IF (IGET(214)>0) THEN
          IF (LVLS(LP,IGET(214))>0) THEN
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=QS1(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(214))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(214))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
         ENDIF
!
!---  GRAUPEL
         IF (IGET(255)>0) THEN
          IF (LVLS(LP,IGET(255))>0) THEN
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=QG1(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(255))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(255))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
         ENDIF
!
!---  TOTAL CONDENSATE
         IF (IGET(215)>0) THEN
          IF (LVLS(LP,IGET(215))>0) THEN 
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=C1D(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(215))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(215))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
         ENDIF
!
!    TOTAL CLOUD COVER
         IF (IGET(222)>0) THEN
          IF (LVLS(LP,IGET(222))>0) THEN 
             DO J=JSTA,JEND
             DO I=ISTA,IEND
               GRID1(I,J)=CFRSIG(I,J)
             ENDDO
             ENDDO
	  if(grib=="grib2" )then
            cfld=cfld+1
            fld_info(cfld)%ifld=IAVBLFLD(IGET(222))
            fld_info(cfld)%lvl=LVLSXML(LP,IGET(222))
            datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=GRID1(ista:iend,jsta:jend)
          endif
          ENDIF
         ENDIF
!***  END OF MAIN VERTICAL LOOP
!     
  310   CONTINUE
!***  ENDIF FOR IF TEST SEEING IF WE WANT ANY OTHER VARIABLES
!
      ENDIF
!
!
!     
!     END OF ROUTINE.
!
      RETURN
      END