    SUBROUTINE SLPSIGSPLINE(PD,FIS,TSIG,QSIG &
    ,            SPL,LSL,DETA,PT,PSLP)

!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
!   SUBROUTINE:  SLPSIGSPLINE      SLP REDUCTION
!   02-1981 TO 03 1998 JANJIC, D. JOVIC, S. NICKOVIC
!   09-2000  H CHUANG: MODIFIED TO BE INCLUDED IN OPERATIONAL
!                      MPI QUILT
!   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

! ABSTRACT:  THIS ROUTINE COMPUTES THE SEA LEVEL PRESSURE
!            REDUCTION USING CUBIC SPLINE FITTING
!            METHOD OR THE STANDARD NCEP REDUCTION FOR
!            SIGMA COORDINATES.

! USAGE:  CALL SLPSIG FROM SUBROUITNE QUILT

!   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

!   SUBPROGRAMS CALLED:
!     UNIQUE:
!             NONE

!-----------------------------------------------------------------------
    INCLUDE "parmeta.f90"
    INCLUDE "COMM_PARA.f90"
    INCLUDE "mpif.h"
    INCLUDE "mpp.h"
    PARAMETER (LMP1=LM+1)
! p      parameter(nsmud=100,lp2=lm+2,lsmp2=lsm+2)
    parameter(nsmud=25,lp2=lm+2,lsmp2=lsm+2)
    parameter(iii=22,jjj=293)
!-----------------------------------------------------------------------
    PARAMETER &
    (RD=287.04,ROG=RD/9.8,PQ0=379.90516,A2=17.2693882 &
    , A3=273.16,A4=35.86,GAMMA=6.5E-3,RGAMOG=GAMMA*ROG &
    , H1M12=1.E-12,g=9.8)
!-----------------------------------------------------------------------
    REAL :: &
    PD(IM,MY_JSD:MY_JED),FIS(IM,MY_JSD:MY_JED) &
    ,tmask(IM,MY_JSD:MY_JED) &
    ,FI(IM,MY_JSD:MY_JED,LSMP2),TSIG(IM,MY_JSD:MY_JED,LM) &
    ,QSIG(IM,MY_JSD:MY_JED,LM)
    REAL :: &
    PSLP(IM,MY_JSD:MY_JED),fsll(IM,MY_JSD:MY_JED)
    REAL :: &
    SPL(LSM)
    REAL :: &
    DETA(LM),AETA(LM),ETA(LMP1),zth(lm+2),hcol  (lm+2) &
    ,y2(lm+2),ztsl  (lsmp2),ovrlx (lsmp2),hcolsl(lsmp2) &
    ,phld  (lm+2),qhld  (lm+2) &
    ,hcol3(3),pcol3(3),y23(3)
!-----------------------------------------------------------------------
    data &
    y2/lp2*0./
!-----------------------------------------------------------------------
    data &
    ovrlx/lsmp2*0.175/
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
    INTEGER :: &
    IHE(JM),IHW(JM)
!-----------------------------------------------------------------------
    LOGICAL :: &
    SIGMA,STDRD
!-----------------------------------------------------------------------
    STDRD= .FALSE. 
!-----------------------------------------------------------------------
!***
!***  CALCULATE THE I-INDEX EAST-WEST INCREMENTS
!***
!      print*,'in slpsig'
    DO J=1,JM
        IHE(J)=MOD(J+1,2)
        IHW(J)=IHE(J)-1
    ENDDO
!-----------------------------------------------------------------------
!***
!***  INITIALIZE ARRAYS.  LOAD SLP ARRAY WITH SURFACE PRESSURE.
!***
!C!$omp parallel do
    DO J=JSTA_I,JEND_I
        DO I=1,IM
            PSLP(I,J)=PD(I,J)+PT
        ENDDO
    ENDDO

    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))
    ENDDO

    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)
            endif
            y2(lmd)=0.
        !	 if(i.eq.iii.and.j.eq.jjj)then
        !	  print*,'sample input height'
        !	  do l=1,lmd
        !	   print*,l,exp(zth(l)**0.5),hcol(l)
        !	  end do
        !	 end if
            call splinef(lmd,zth,hcol,y2,lsmp2,ztsl,hcolsl,phld,qhld)
            do 240 l=1,lsmp2
                fi(i,j,l)=hcolsl(l)
            !	  if(i.eq.iii.and.j.eq.jjj)print*,'sample H'
            !     +    ,exp(ZTSL(L)**0.5),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.
                endif
                if(href >= fi(i,j,l))then
                    tmask(i,j)=rlx
                else
                    tmask(i,j)=0.
                endif
            311 END DO
        310 END DO
    !        print*,'before relaxation at',l
        do 320 n=1,nsmud
        !	 if(l.eq.41)print*,'calling update at n=',n
            call update(fi(1,my_jsd,l))
        !         if(l.eq.41)print*,'after calling update at n=',n
        !	 print*,'jsta_im jend_im=',jsta_im,jend_im
            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.
                    endif
                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)
                    endif
                341 END DO
            340 END DO
        320 END DO
    !        if(i.eq.iii.and.j.eq.jjj)print*,'H after relaxation',i,j
    !     +,  fi(i,j,l)
    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
                    gar=exp(ztsl(lsmp2)**0.5)+5000.     !assume pslp < 1100 mb
                    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
                            go to 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
            endif
        !         if(i.eq.iii.and.j.eq.jjj)then
        !	  print*,'on P after relaxation at ',iii,jjj
        !	  do l=1,lsmp2
        !	   print*,'l,H,P= ',l,fi(i,j,l),exp(ztsl(l)**0.5)
        !          end do
        !	  print*,'pslp = ',pslp(i,j)
        ! 	 end if
        !         if(pslp(i,j).gt.100200. .and. pslp(i,j).lt.100400.)
        !     +   print*,'target pslp',i,j,pslp(i,j)
        601 END DO
    600 END DO
!--------------------------------------------------------------------
!     SKIP THE STANDARD SCHEME.
!--------------------------------------------------------------------
    GO TO 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
                ENDIF
                PSLP(I,J)=EXP(SLPP)
            ENDIF
    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)

!     ******************************************************************
!     *                                                                *
!     *  THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE        *
!     *  PROGRAMED FOR A SMALL SCALAR MACHINE.                         *
!     *                                                                *
!     *  PROGRAMER[ Z. JANJIC, YUGOSLAV FED. HYDROMET. INST., BEOGRAD  *
!     *                                                                *
!     *                                                                *
!     *                                                                *
!     *  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.                *
!     *                                                                *
!     ******************************************************************

    DIMENSION &
    XOLD(NOLD),YOLD(NOLD),Y2(NOLD),P(NOLD),Q(NOLD) &
    ,XNEW(NNEW),YNEW(NNEW)
!-----------------------------------------------------------------------
    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) GO TO 700
!-----------------------------------------------------------------------
    K=3

    100 DXL=DXR
    DYDXL=DYDXR
    DXR=XOLD(K+1)-XOLD(K)
!      if(i.eq.120.and.j.eq.279)then
!       print*,'in spline',k,XOLD(K),XOLD(K+1),dxr
!      end if
    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) GO TO 100
!-----------------------------------------------------------------------
    700 K=NOLDM1

    200 Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)

    K=K-1
    IF(K > 1) GO TO 200
!-----------------------------------------------------------------------
    K1=1

    300 XK=XNEW(K1)

    DO 400 K2=2,NOLD
        IF(XOLD(K2) <= XK) GO TO 400
        KOLD=K2-1
        GO TO 450
    400 END DO
    YNEW(K1)=YOLD(NOLD)
    GO TO 600

    450 IF(K1 == 1)   GO TO 500
    IF(K == KOLD) GO TO 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) GO TO 300
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE SPLINEF
