    SUBROUTINE SLPSIGSPLINE(PD, FIS, TSIG, QSIG, SPL, LSL, DETA, PT, PSLP)
!>--------------------------------------------------------------------------------------------------
!> SUBROUTINE SLPSIGSPLINE
!>
!> SUBROUTINE: SLPSIGSPLINE - SLP REDUCTION
!> PROGRAMMER: JANJIC, D. JOVIC, S. NICKOVIC
!> ORG: ?????
!> DATE: 81-02-??
!>
!> PROGRAMMER: JANJIC, D. JOVIC, S. NICKOVIC
!> ORG: ?????
!> DATE: 81-02-??
!>
!> ABSTRACT: THIS ROUTINE COMPUTES THE SEA LEVEL PRESSURE REDUCTION USING CUBIC SPLINE FITTING
!>           METHOD OR THE STANDARD NCEP REDUCTION FOR SIGMA COORDINATES.
!>
!> PROGRAM HISTORY LOG:
!> 81-02-??  JANJIC       - ORIGINATOR
!> 00-09-??  H CHUANG     - INCLUDED IN OPERATIONAL MPI QUILT 
!> 18-01-15  LUCCI        - MODERNIZATION OF THE CODE, INCLUDING:
!>                          * F77 TO F90/F95
!>                          * INDENTATION & UNIFORMIZATION CODE
!>                          * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!>                          * DOCUMENTATION WITH DOXYGEN
!>                          * OPENMP FUNCTIONALITY
!>
!> NOTE: CURRENT POST ONLY PROCESS UP TO 1000 MB PRESSURE LEVEL 2 ADDITIONAL LEVELS (1025 1050) ARE
!>       ADDED DURING SPLINE FITTING COMPUTATION BUT THE FIELDS ON THESE TWO LEVELS ARE NOT OUTPUT 
!>       TO GRID FILES.
!>
!> INPUT ARGUMENT LIST:
!> PD   - SFC PRESSURE MINUS PTOP
!> FIS  - SURFACE GEOPOTENTIAL
!> TSIG - TEMPERATURE
!> QSIG - SPECIFIC HUMIDITY
!> FI   - GEOPOTENTIAL
!> PT   - TOP PRESSURE OF DOMAIN
!>
!> OUTPUT ARGUMENT LIST:
!> PSLP - THE FINAL REDUCED SEA LEVEL PRESSURE ARRAY
!>
!> INPUT/OUTPUT ARGUMENT LIST:
!> NONE
!>
!> USE MODULES: F77KINDS
!>              GLB_TABLE
!>              MAPPINGS
!>              MPPCOM
!>              PARA
!>              PARMETA
!>              TEMPCOM
!>              TOPO
!>
!> DRIVER     : QUILT
!>
!> CALLS      : SPLINEF
!>              UPDATE 
!>--------------------------------------------------------------------------------------------------
    USE F77KINDS
    USE GLB_TABLE
    USE MAPPINGS
    USE MPPCOM
    USE PARA
    USE PARMETA
    USE TEMPCOM
    USE TOPO
!
    IMPLICIT NONE	
!
    INCLUDE "mpif.h"
!
    INTEGER(KIND=I4KIND), PARAMETER :: LMP1   = LM  + 1
    INTEGER(KIND=I4KIND), PARAMETER :: NSMUD  =  25 
    INTEGER(KIND=I4KIND), PARAMETER :: LP2    = LM  + 2
    INTEGER(KIND=I4KIND), PARAMETER :: LSMP2  = LSM + 2
    INTEGER(KIND=I4KIND), PARAMETER :: III    =  22
    INTEGER(KIND=I4KIND), PARAMETER :: JJJ    = 293
!
    REAL   (KIND=R4KIND), PARAMETER :: RD     = 287.04
    REAL   (KIND=R4KIND), PARAMETER :: ROG    = RD / 9.8
    REAL   (KIND=R4KIND), PARAMETER :: PQ0    = 379.90516
    REAL   (KIND=R4KIND), PARAMETER :: A2     =  17.2693882 
    REAL   (KIND=R4KIND), PARAMETER :: A3     = 273.16
    REAL   (KIND=R4KIND), PARAMETER :: A4     =  35.86
    REAL   (KIND=R4KIND), PARAMETER :: GAMMA  =   6.5E-3
    REAL   (KIND=R4KIND), PARAMETER :: RGAMOG = GAMMA * ROG
    REAL   (KIND=R4KIND), PARAMETER :: H1M12  =   1.E-12
    REAL   (KIND=R4KIND), PARAMETER :: G      =   9.8
!
    REAL   (KIND=R4KIND), DIMENSION(IM, MY_JSD:MY_JED)                    , INTENT(IN)          ::&
    & PD      , FIS     
!
    REAL   (KIND=R4KIND), DIMENSION(IM, MY_JSD:MY_JED)                                          ::&
    & TMASK   , FSLL
!
    REAL   (KIND=R4KIND), DIMENSION(IM, MY_JSD:MY_JED, LSMP2)                                   ::&
    & FI
!
    REAL   (KIND=R4KIND), DIMENSION(IM, MY_JSD:MY_JED, LM)                , INTENT(IN)          ::&
    & TSIG    , QSIG
!
    REAL   (KIND=R4KIND), DIMENSION(IM, MY_JSD:MY_JED)                    , INTENT(OUT)         ::&
    & PSLP
!
    REAL   (KIND=R4KIND), DIMENSION(LSM)                                  , INTENT(IN)          ::&
    & SPL
!
    REAL   (KIND=R4KIND), DIMENSION(LM)                                   , INTENT(IN)          ::&
    & DETA 
!
    REAL   (KIND=R4KIND), DIMENSION(LM)                                                         ::&
    & AETA
!
    REAL   (KIND=R4KIND), DIMENSION(LMP1)                                                       ::&
    & ETA
!
    REAL   (KIND=R4KIND), DIMENSION(LM+2)                                                       ::&
    & ZTH     , HCOL    , Y2      , PHLD    , QHLD
!
    REAL   (KIND=R4KIND), DIMENSION(LSMP2)                                                      ::&
    & ZTSL    , OVRLX   , HCOLSL
!
    REAL   (KIND=R4KIND), DIMENSION(3)                                                          ::&
    & HCOL3   , PCOL3   , Y23
!
    DATA Y2 /LP2 * 0./
    DATA OVRLX /LSMP2 * 0.175/
!
    INTEGER(KIND=I4KIND), DIMENSION(JM)                                                         ::&
    & IHE     , IHW
!
    LOGICAL(KIND=L4KIND)                                                                        ::&
    & SIGMA, STDRD
!------------------------
! IMPLICIT NONE VARIABLES
!------------------------
    INTEGER(KIND=I4KIND)                                                  , INTENT(IN)          ::&
    & LSL
!
    REAL   (KIND=R4KIND)                                                  , INTENT(IN)          ::&
    & PT
!
    INTEGER(KIND=I4KIND)                                                                        ::&
    & I       , J       , L       , IHL     , IHH     , IVI     , LMD     , N       , LLL     ,   &
    & LMA
!
    REAL   (KIND=R4KIND)                                                                        ::&
    & ZTBOT   , PDP     , ALP1L   , ALP2L   , DH      , X       , DUM     , D2      , RLX     ,   &
    & HREF    , GAR     , SLP1    , ALP1U   , ALP2U   , ALPP1   , SLOP    , SLPP    , TTT
!
    STDRD = .FALSE. 
!------------------------------------------ 
!CALCULATE THE I-INDEX EAST-WEST INCREMENTS
!------------------------------------------ 
    DO J=1,JM
        IHE(J) = MOD(J+1,2)
        IHW(J) = IHE(J)-1
    END DO
!---------------------------------------------------------- 
! INITIALIZE ARRAYS.  LOAD SLP ARRAY WITH SURFACE PRESSURE.
!----------------------------------------------------------
!
!$omp parallel do
!
    DO J=JSTA_I,JEND_I
        DO I=1,IM
            PSLP(I,J) = PD(I,J) + PT
        END DO
    END DO
!
    DO 100 L=1,LSM
        ZTSL(L) = ALOG(SPL(L) ) ** 2
100 END DO
!
    ZTSL(LSM+1) = ALOG(102500.) ** 2
    ZTSL(LSMP2) = ALOG(105000.) ** 2
!
    ZTBOT = ZTSL(LSMP2)
!---------------------  
! COMPUTE ETA AND AETA
!--------------------- 
    ETA(1) = 0.0
!
    DO L=2,LM+1
         ETA(L) =        ETA(L-1) + DETA(L-1)
    END DO
!
    DO L=1,LM
        AETA(L) = 0.5 * (ETA(L  ) +  ETA(L+1))
    END DO
!
    DO 200 J=JSTA_I,JEND_I
        IHL = 1
        IHH = IM-1 + MOD(J,2)
!
        DO 201 I=IHL,IHH
            PDP        =  PD(I,J)
            HCOL(LM+1) = FIS(I,J) / G
!
            ALP1L      =  ALOG(PT  + PDP)
            ZTH(LM+1)  = (ALOG(PDP + PT )) ** 2
!
            DO 220 IVI=1,LM
                L       = LM + 1 - IVI
                ALP1U   =  ALOG(ETA(L) * PDP + PT)
                ALP2U   = (ALOG(ETA(L) * PDP + PT)) ** 2
                DH      = (QSIG(I,J,L) * 0.608 + 1.) * (ALP1L - ALP1U) * TSIG(I,J,L) * ROG
!
                 ZTH(L) = ALP2U
                HCOL(L) = HCOL(L+1) + DH
!
                ALP1L   = ALP1U
        220 END DO
!
            IF (ZTH(LM+1) >= ZTBOT) THEN
                LMD = LM + 1
            ELSE
                LMD        = LM + 2
                ZTH(LM+2)  = ZTBOT
                X          = ZTBOT - ZTH(LM+1)
                DUM        = (HCOL(LM+1) - HCOL(LM-5)) / (ZTH(LM+1) - ZTH(LM-5))
                D2         = 0.
                HCOL(LM+2) = D2 * X * X + DUM * X + HCOL(LM+1)
            END IF
!
            Y2(LMD) = 0.
!
            CALL SPLINEF(LMD, ZTH, HCOL, Y2, LSMP2, ZTSL, HCOLSL, PHLD, QHLD)
!
            DO 240 L=1,LSMP2
                FI(I,J,L) = HCOLSL(L)
        240 END DO
!
    201 END DO
!
 200 END DO
!-------------------------------------------------- 
! FILLING REMAINING UNDEFINED MASS POINT BOUNDARIES 
!-------------------------------------------------- 
    DO 300 L=1,LSMP2
!
        IF (MOD(JSTA_I,2) < 1) THEN !EVEN JSTA_I
!
            DO 351 J=JSTA_I,JEND_I,2
                FI(IM,J,L) = FI(IM-1,J,L)
        351 END DO
!
        ELSE
!
            DO 352 J=JSTA_I+1,JEND_I,2
                FI(IM,J,L) = FI(IM-1,J,L)
            352 END DO
!
        END IF
!
        RLX = OVRLX(L)
!
        DO 310 J=JSTA_I,JEND_I
            IHL = 1
            IHH = IM - 1 + MOD(J,2)
!
            DO 311 I=IHL,IHH
                HREF = FIS(I,J) / G
!
                IF (HREF > 300.) THEN
                    HREF = HREF + 1500.
                END IF
!
                IF (HREF >= FI(I,J,L)) THEN
                    TMASK(I,J) = RLX
                ELSE
                    TMASK(I,J) = 0.
                END IF
!
        311 END DO
!
    310 END DO
!
        DO 320 N=1,NSMUD
!
            CALL UPDATE(FI(1,MY_JSD,L))
!
            DO 330 J=JSTA_IM,JEND_IM
                IHL = 1 + MOD(J,2)
                IHH = IM-1
!
                DO 331 I=IHL,IHH
                    IF (TMASK(I,J) > 0.05) THEN
                        FSLL(I,J) = FI(I + IHW(J),J-1,L) + FI(I + IHE(J),J-1,L)                   &
    &                             + FI(I + IHW(J),J+1,L) + FI(I + IHE(J),J+1,L) - FI(I,J,L) * 4.
                    END IF
            331 END DO
!
        330 END DO
!------------------------- 
! PRINT*,'AFTER SMOOTHING'
!------------------------- 
            DO 340 J=JSTA_IM,JEND_IM
                IHL = 1 + MOD(J,2)
                IHH = IM-1
!
                DO 341 I=IHL,IHH
                    IF (TMASK(I,J) > 0.05) THEN
                        FI(I,J,L) = FSLL(I,J) * TMASK(I,J) + FI(I,J,L)
                    END IF
            341 END DO
!
        340 END DO
!
    320 END DO
!
300 END DO
!-------------------  
! SEA LEVEL PRESSURE 
!-------------------
    DO 600 J=JSTA_I,JEND_I
        IHL = 1
        IHH = IM - 1 + MOD(J,2)
!
        DO 601 I=IHL,IHH
            IF (FIS(I,J) > -1. .AND. FIS(I,J) < 1.) THEN
                PSLP(I,J) = PT + PD(I,J)
            ELSE
                IF (FI(I,J,LSMP2) > 0) THEN
!---------------------- 
! ASSUME PSLP < 1100 MB
!---------------------- 
                    GAR = EXP(ZTSL(LSMP2) ** 0.5) + 5000. 
                    HCOL3(1) = FI(I,J,LSMP2) + (FI(I,J,LSMP2) - FI(I,J,LSMP2-5))                  &
    &                        * (ALOG(GAR) ** 2. - ZTSL(LSMP2)) / (ZTSL(LSMP2) - ZTSL(LSMP2-5))
                    HCOL3(2) = FI(I,J,LSMP2  )
                    HCOL3(3) = FI(I,J,LSMP2-1)
                    PCOL3(1) = GAR
                    PCOL3(2) = EXP(ZTSL(LSMP2  ) ** 0.5)
                    PCOL3(3) = EXP(ZTSL(LSMP2-1) ** 0.5)
                ELSE
                    DO L=LSMP2,2,-1
                        IF (FI(I,J,L) <= 0. .AND. FI(I,J,L-1) > 0.) THEN
                            LLL = L
                            GOTO 603
                        END IF
                    END DO
!
                603 HCOL3(1) = FI(I,J,LLL  )
                    HCOL3(2) = FI(I,J,LLL-1)
                    HCOL3(3) = FI(I,J,LLL-2)
!
                    PCOL3(1) = EXP(ZTSL(LLL  ) ** 0.5)
                    PCOL3(2) = EXP(ZTSL(LLL-1) ** 0.5)
                    PCOL3(3) = EXP(ZTSL(LLL-2) ** 0.5)
                END IF
!
                DO L=1,3
                    Y23(L) = 0.0
                END DO
!
                IF (HCOL3(1) == HCOL3(2) .OR. HCOL3(2) == HCOL3(3))                               &
    &               PRINT*, 'PROBLEM SLP', I, J, HCOL3(1), HCOL3(2), HCOL3(3)
!
                CALL SPLINEF(3, HCOL3, PCOL3, Y23, 1, 0.0, SLP1, PHLD, QHLD)
                PSLP(I,J) = SLP1
            END IF
!
    601 END DO
!
600 END DO
!-------------------------- 
! SKIP THE STANDARD SCHEME.
!-------------------------- 
    GOTO 430
!-------------------------------------------------------------------------
! IF YOU WANT THE "STANDARD" ETA/SIGMA REDUCTION THIS IS WHERE IT IS DONE.
!-------------------------------------------------------------------------
    DO 410 J=JSTA_I,JEND_I
        DO 410 I=1,IM
            IF (FIS(I,J) >= 1.) THEN
                LMA   = LM
                ALPP1 = ALOG(PD(I,J) + PT)
                SLOP  = 0.0065 * ROG * TSIG(I,J,LM)
!
                IF (SLOP < 0.50) THEN
                    SLPP = ALPP1 + FIS(I,J) / (RD * TSIG(I,J,LMA))
                ELSE
                    TTT = -(ALOG(PD(I,J) + PT) + ALPP1) *   SLOP * 0.50 + TSIG(I,J,LMA)
                    SLPP = (-TTT + SQRT(TTT * TTT + 2. * SLOP * (FIS(I,J) / RD                    &
    &                    +  (TTT + 0.50 * SLOP * ALPP1) * ALPP1))) / SLOP
                END IF
!
                PSLP(I,J) = EXP(SLPP)
!
            END IF
410 END DO
!--------------------------------------------------------------------------------------------------
! AT THIS POINT WE HAVE A SEA LEVEL PRESSURE FIELD BY EITHER METHOD. 5-POINT AVERAGE THE FIELD ON
! THE E-GRID.
!--------------------------------------------------------------------------------------------------
!
430 CONTINUE
!
    RETURN
!
    END SUBROUTINE SLPSIGSPLINE
!
!
!
    SUBROUTINE SPLINEF(NOLD, XOLD, YOLD, Y2, NNEW, XNEW, YNEW, P, Q)
!>--------------------------------------------------------------------------------------------------
!> SUBROUTINE SPLINEF
!>
!> SUBROUTINE: SLPSIGSPLINE - SLP REDUCTION
!> PROGRAMMER: JZ. JANJIC, YUGOSLAV FED. HYDROMET.
!> ORG: INST. BEOGRAD
!> DATE: 81-02-??
!>
!> ABSTRACT:
!> THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE PROGRAMED FOR A SMALL SCALAR MACHINE.
!>
!> PROGRAM HISTORY LOG:
!> 81-02-??  JANJIC       - ORIGINATOR
!> 18-01-15  LUCCI        - MODERNIZATION OF THE CODE, INCLUDING:
!>                          * F77 TO F90/F95
!>                          * INDENTATION & UNIFORMIZATION CODE
!>                          * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!>                          * DOCUMENTATION WITH DOXYGEN
!>                          * OPENMP FUNCTIONALITY
!>
!> INPUT ARGUMENT LIST:
!> NOLD - 
!> XOLD -
!> YOLD -
!> NNEW -
!> XNEW - 
!> YNEW -
!>
!> OUTPUT ARGUMENT LIST:
!> NONE
!>
!> INPUT/OUTPUT ARGUMENT LIST:
!> Y2   -
!> P    -
!> Q    - 
!>
!> USE MODULES: F77KINDS
!>
!> DRIVER     : SLPSIGSPLINE
!>
!> CALLS      : -----
!--------------------------------------------------------------------------------------------------
!
!--------------------------------------------------------------------------------------------------
! NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3.
! XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE FUNCTION ARE GIVEN. MUST BE IN 
!        ASCENDING ORDER.
! YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. 
! Y2   - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL SPLINE IS FITTED Y2(1)=0. 
!        AND Y2(NOLD)=0. MUST BE SPECIFIED.
! NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED.
! XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE FUNCTION ARE CALCULATED. XNEW(K) 
!        MUST BE GE XOLD(1) AND LE XOLD(NOLD).
! YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2
!--------------------------------------------------------------------------------------------------
    USE F77KINDS
!
    IMPLICIT NONE
!
    INTEGER(KIND=I4KIND)                                                  , INTENT(IN)          ::&
    & NOLD    , NNEW
!
    REAL   (KIND=R4KIND), DIMENSION(NOLD)                                 , INTENT(IN)          ::&
    & XOLD    , YOLD    
!
    REAL   (KIND=R4KIND), DIMENSION(NOLD)                                 , INTENT(INOUT)       ::&
    & Y2      , P       , Q
!
    REAL   (KIND=R4KIND), DIMENSION(NNEW)                                 , INTENT(IN)          ::&
    & XNEW     
!
    REAL   (KIND=R4KIND), DIMENSION(NNEW)                                 , INTENT(INOUT)       ::&
    & YNEW
!------------------------
! IMPLICIT NONE VARIABLES
!------------------------
    INTEGER(KIND=I4KIND)                                                                        ::&
    & NOLDM1  , K       , K1      , K2      , KOLD     
!
    REAL   (KIND=R4KIND)                                                                        ::&
    & DXL     , DXR     , DYDXL   , DYDXR   , RTDXC   , DXC     , DEN     , XK      , Y2K     ,   &
    & Y2KP1   , DX      , RDX     , AK      , BK      , CK      , X       , XSQ
!
    NOLDM1 = NOLD - 1
!
    DXL   =  XOLD(2)  - XOLD(1)
    DXR   =  XOLD(3)  - XOLD(2)
    DYDXL = (YOLD(2)  - YOLD(1)) / DXL
    DYDXR = (YOLD(3)  - YOLD(2)) / DXR
    RTDXC = .5 / (DXL + DXR)
!
    P(1) =  RTDXC * (6. * (DYDXR - DYDXL) - DXL * Y2(1))
    Q(1) = -RTDXC * DXR
!
    IF (NOLD == 3) GOTO 700
!
    K = 3
!
100 DXL   = DXR
    DYDXL = DYDXR
    DXR   = XOLD(K+1) - XOLD(K)
!
    DYDXR = (YOLD(K+1) - YOLD(K)) / DXR
    DXC   = DXL + DXR
    DEN   = 1. / (DXL * Q(K-2) + DXC + DXC)
!
    P(K-1) =  DEN * (6. * (DYDXR - DYDXL) - DXL * P(K-2))
    Q(K-1) = -DEN * DXR
!
    K = K + 1
!
    IF (K < NOLD) GOTO 100
!
700 K = NOLDM1
!
200 Y2(K) = P(K-1) + Q(K-1) * Y2(K+1)
!
    K = K - 1
!
    IF (K > 1) GOTO 200
!
    K1 = 1
!
300 XK = XNEW(K1)
!
    DO 400 K2=2,NOLD
        IF (XOLD(K2) <= XK) GOTO 400
        KOLD = K2 - 1
        GOTO 450
400 END DO
!
    YNEW(K1) = YOLD(NOLD)
!
    GOTO 600
!
450 IF (K1 == 1)   GOTO 500
!
    IF (K == KOLD) GOTO 550
!
500 K = KOLD
!
    Y2K   = Y2(K  )
    Y2KP1 = Y2(K+1)
!
    DX    = XOLD(K+1) - XOLD(K)
    RDX   = 1. / DX
!
    AK = .1666667 * RDX * (Y2KP1 - Y2K)
    BK = .5 * Y2K
    CK = RDX * (YOLD(K+1) - YOLD(K)) - .1666667 * DX * (Y2KP1 + Y2K + Y2K)
!
550 X   = XK - XOLD(K)
    XSQ = X * X
!
    YNEW(K1) = AK * XSQ * X + BK * XSQ + CK * X + YOLD(K)
!
600 K1 = K1 + 1
!
    IF (K1 <= NNEW) GOTO 300
!
    RETURN
!
    END SUBROUTINE SPLINEF
