       subroutine slp

       implicit none
       include 'param_o.h'
       include 'vrbls.comm'
       include 'extra.comm'
       include 'mapot.comm'
       include 'const.h'
       include 'dynam.comm'
       include 'masks.comm'
       include 'params'

       integer,parameter::nrlx1=500
       integer::i,j,n,l,lhmnt,lmap1,ll,nrlx,kmn,kmntm,lmst,li,kmm, &
                km,n1
       integer::isCorner2,lma
       real,dimension(0:im+1,0:jm+1,nm)::ttv,pdsl1,pbi,slpx
       integer,dimension((im+2)*(jm+2)*nm)::imnt,jmnt,nmnt
       real::tgss,pbin,phbi,ptin,trtv,phti,dposp,prbin,prtin,alpp1,slop,  &
             slpp,ttt
       real,parameter::OVERRC=1.50,AD05=OVERRC*0.05,CFT0=OVERRC-1., &
                       r=287.04,rog=r/9.8
       logical:: stdrd

       if(sigma)then
         stdrd=.true.
       else
         stdrd=.false.
       endif
!
!-----------------------------------------------------------------------
!
!***  FIND THE MOST ELEVATED GLOBAL LAYER WHERE THERE IS LAND
!
!-----------------------------------------------------------------------
        DO L=1,LM
          do n=1,nm
          DO J=1,jm
          DO I=1,IM
            IF(htm(I,J,n,l).LT.0.5)GO TO 666
          ENDDO
          ENDDO
          ENDDO
!
!***  IF WE GET TO HERE, WE HAVE ALL ATMOSPHERE
!
          GOTO 667
  666     CONTINUE
!
!***  IF WE GET HERE, WE HAVE FOUND A NON_ATM POINT
!
          LHMNT=L
          GOTO 668
  667     CONTINUE
        ENDDO
!
        LHMNT=LM+1
!cccccc GO TO 669
  668   CONTINUE

!-----------------------------------------------------------------------
!***
!***  INITIALIZE ARRAYS.  LOAD SLP ARRAY WITH SURFACE PRESSURE.
!***
      do n=1,nm
      DO J=1,jm
      DO I=1,IM
        PSLP(I,J,n)=0.
        TTV(I,J,n)=0.
      ENDDO
      ENDDO
      ENDDO
!
      do 110 n=1,nm
      DO 110 J=1,jm
      DO 110 I=1,IM
      PDSL1(I,J,n)=RES(I,J,n)*PD(I,J,n)
      PSLP(I,J,n)=PD(I,J,n)+PT
      PBI (I,J,n)=PSLP(I,J,n)
  110 CONTINUE

      if(.not.stdrd)then

      LL=LM
!
      do n=1,nm
      DO J=1,jm
      DO I=1,im
        IF(HTM(I,J,n,LL).LE.0.5) THEN
          TGSS=FIS(I,J,n)/(R*ALOG((PDSL1(I,J,n)+PT)/(PD(I,J,n)+PT)))
          LMAP1=LMH(I,J,n)+1
!
	  if(n.eq.3.and.i.eq.30.and.j.eq.101)then
!	    print *,i,j,tgss,fis(i,j,n),r,pdsl1(i,j,n),pd(i,j,n)
!	    stop
          endif
          DO 260 L=LMAP1,LM
          T(I,J,n,L)=TGSS
  260     CONTINUE
!
        ENDIF
      ENDDO
      ENDDO
      ENDDO

!----------------------------------------------------------------
!
!***  CREATE A TEMPORARY TV ARRAY, AND FOLLOW BY SEQUENTIAL
!***  OVERRELAXATION, DOING NRLX PASSES.
!
!----------------------------------------------------------------
      NRLX=NRLX1
!
!----------------------------------------------------------------
!----------------------------------------------------------------
      DO 300 L=LHMNT,LM
!----------------------------------------------------------------
!----------------------------------------------------------------
!
      KMN=0
      KMNTM=0
!
      do 240 n=1,nm
      DO 240 J=1,jm
      DO 240 I=1,im
      IF(HTM(I,J,n,L).GT.0.5)GO TO 240
      KMN=KMN+1
      IMNT(KMN)=I
      JMNT(KMN)=J
      nmnt(kmn)=n
  240 CONTINUE
!
      KMNTM=KMN
!
      do 270 n=1,nm
      DO 270 J=1,jm
      DO 270 I=1,IM
      TTV(I,J,n)=T(I,J,n,L)
  270 CONTINUE
!
!----------------------------------------------------------------
!***  FOR GRID BOXES NEXT TO MOUNTAINS REPLACE TTV BY AN "EQUIVALENT"
!***  TV, ONE WHICH CORRESPONDS TO THE CHANGE IN P BETWEEN REFERENCE
!***  INTERFACE GEOPOTENTIALS, INSTEAD OF BETWEEN LAYER INTERFACES
!----------------------------------------------------------------
!
      do n=1,nm
      DO J=1,jm
      DO I=1,im
        IF(HTM(I,J,n,L).GT.0.5.AND.  &
           HTM(I,J-1,n,L)*HTM(I+1,J-1,n,L)  &
          *HTM(I,J+1,n,L)*HTM(I+1,J+1,n,L)  &
          *HTM(I-1     ,J  ,n,L)*HTM(I+1     ,J  ,n,L)  &
          *HTM(I-1     ,J-1,n,L)*HTM(I-1     ,J+1,n,L).LT.0.5)THEN
          LMST=LMH(I,J,n)
!***
!***  FIND P AT THE REFERENCE INTERFACE GEOPOTENTIAL AT THE BOTTOM
!***
          PBIN=PT+PD(I,J,n)
          PHBI=DFL(LMST+1)
!
          DO LI=LMST,1,-1
            PTIN=PBIN-DETA(LI)*PD(I,J,n)*RES(I,J,n)
            TRTV=2.*R*T(I,J,n,LI)*(1.+0.608*Q(I,J,n,LI))
            PHTI=PHBI+TRTV*(PBIN-PTIN)/(PBIN+PTIN)
	  if(i.eq.17.and.j.eq.89.and.n.eq.5)then
!	    print *,l,pbin,ptin,deta(li),res(i,j,n),pd(i,j,n)
          endif
            IF(PHTI.GE.DFL(L+1))GO TO 273
            PBIN=PTIN
            PHBI=PHTI
          ENDDO
!
  273     DPOSP=(PHTI-DFL(L+1))/TRTV
          PRBIN=(1.+DPOSP)/(1.-DPOSP)*PTIN

!***
!***  FIND P AT THE REFERENCE INTERFACE GEOPOTENTIAL AT THE TOP
!***
          PBIN=PT+PD(I,J,n)
          PHBI=DFL(LMST+1)
!
          DO LI=LMST,1,-1
            PTIN=PBIN-DETA(LI)*PD(I,J,n)*RES(I,J,n)
            TRTV=2.*R*T(I,J,n,LI)*(1.+0.608*Q(I,J,n,LI))
            PHTI=PHBI+TRTV*(PBIN-PTIN)/(PBIN+PTIN)
            IF(PHTI.GE.DFL(L))GO TO 275
            PBIN=PTIN
            PHBI=PHTI
          ENDDO
!
  275     DPOSP=(PHTI-DFL(L))/TRTV
          PRTIN=(1.+DPOSP)/(1.-DPOSP)*PTIN
!
          TTV(I,J,n)=(DFL(L)-DFL(L+1))/(2.*R)*(PRBIN+PRTIN)/(PRBIN-PRTIN)
        ENDIF
      ENDDO
      ENDDO
      ENDDO

!----------------------------------------------------------------
      KMM=KMNTM
!----------------------------------------------------------------
!***
!***  HERE IS THE RELAXATION LOOP
!***
!----------------------------------------------------------------
      DO 285 N1=1,NRLX
!
      CALL bocoh(TTV,im,jm,nm)   ! Exchange haloes
!
      DO 280 KM=1,KMM
      I=IMNT(KM)
      J=JMNT(KM)
      n=nMNT(KM)
      if(isCorner2(i,j,n,im,jm,nm).eq.0)then
      TTV(I,J,n)=AD05*(4.*(TTV(I,J-1,n)+TTV(I+1,J,n)  &
                        +TTV(I,J+1,n)+TTV(I-1,J,n))  &
                        +TTV(I-1,J-1,n)     +TTV(I-1,J+1,n)  &
                        +TTV(I+1,J-1,n)     +TTV(I+1,J+1,n))  &
                        -CFT0*TTV(I,J,n)
      else
       TTV(i,j,n)=0.25*(TTV(I,J-1,n)+TTV(I+1,J,n)  &
                        +TTV(I,J+1,n)+TTV(I-1,J,n))
      endif
  280 CONTINUE
!
  285 CONTINUE
!----------------------------------------------------------------
!
      if(l.eq.36)then
!        print *,'t(30,101,3),',t(30,101,3,36),ttv(30,101,3)
      endif

      DO 290 KM=1,KMM
      I=IMNT(KM)
      J=JMNT(KM)
      n=nMNT(KM)
      if(l.eq.36)then
!        print *,i,j,n,km
      endif
      T(I,J,n,L)=TTV(I,J,n)
  290 CONTINUE
!
  300 CONTINUE
!     VALUES FOR IMNT AND JMNT ARE FOR LAYER LM - THIS IS WHAT WE WANT
      KMM=KMNTM
!
      DO 320 KM=1,KMM
      I=IMNT(KM)
      J=JMNT(KM)
      n=nMNT(KM)
      LMAP1=LMH(I,J,n)+1
      PBIN=PT+PD(I,J,n)
!
      DO L=LMAP1,LM
        PTIN=PBIN
        DPOSP=(DFL(L)-DFL(L+1))/(2.*R*T(I,J,n,L))
        PBIN=(1.+DPOSP)/(1.-DPOSP)*PTIN
      if((i.eq.35.or.i.eq.36).and.j.eq.58.and.n.eq.3)then
!        print *,i,j,n,l,pd(i,j,n),ptin,dposp
!	stop
      endif

      ENDDO
!
      PSLP(I,J,n)=PBIN
      if(i.eq.70.and.j.eq.186.and.n.eq.4)then
!        print *,i,j,n,lmap1,lm
!	stop
      endif
      if(pslp(i,j,n).gt.123000)then
!        print *,i,j,n,pslp(i,j,n)
!	stop
      endif
  320 CONTINUE
!      stop

      else

        write(6,*) 'doing standard reduction!!!'
      DO 410 n=1,nM
      DO 410 J=1,JM
      DO 410 I=1,IM
      IF(FIS(I,J,n).GE.1.)THEN
        LMA=LMH(I,J,n)
        ALPP1=ALOG(PDSL1(I,J,n)*ETA(LMA+1)+PT)
        SLOP=0.0065*ROG*T(I,J,n,LMA)
        IF(SLOP.LT.0.50)THEN
          SLPP=ALPP1+FIS(I,J,n)/(R*T(I,J,n,LMA))
        ELSE
          TTT=-(ALOG(PDSL1(I,J,n)*ETA(LMA)+PT)+ALPP1)  &
              *SLOP*0.50+T(I,J,n,LMA)
          SLPP=(-TTT+SQRT(TTT*TTT+2.*SLOP*  &
               (FIS(I,J,n)/R+                &
               (TTT+0.50*SLOP*ALPP1)*ALPP1)))/SLOP
        ENDIF
        PSLP(I,J,n)=EXP(SLPP)
      ENDIF
  410 CONTINUE
      endif

      do 440 n=1,nm
      DO 440 J=1,JM
      DO 440 I=1,IM
      SLPX(I,J,n)=PSLP(I,J,n)
  440 CONTINUE
!C
!C
!C!$OMP parallel do private(ihh2)
      do 460 n=1,nm
      DO 460 J=2,JM-1
      DO 460 I=2,im-1
!C
!C***  EXTRA AVERAGING UNDER MOUNTAINS TAKEN OUT, FM, MARCH 96
!C
      SLPX(I,J,n)=0.125*(PSLP(I,J-1,n)+PSLP(I,J+1,n)   &
                     +PSLP(I+1,J,n)+PSLP(I-1,J,n)       &
                     +4.*PSLP(I,J,n))
  460 CONTINUE
!C
!C!$OMP parallel do
      do n=1,nm
      DO J=1,JM
      DO I=1,IM
        PSLP(I,J,n)=SLPX(I,J,n)
      ENDDO
      ENDDO
      ENDDO

      end subroutine slp
