      SUBROUTINE NGMSLP2
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    NGMSLP2     NMC SEA LEVEL PRESSURE REDUCTION
C   PRGRMMR: TREADON         ORG: W/NP2      DATE: 93-02-02       
C     
C ABSTRACT:
C     
C     THIS ROUTINE COMPUTES SEA LEVEL PRESSURE USING THE
C     HYDROSTATIC EQUATION WITH THE SHUELL CORRECTION.  THE
C     FOLLOWING IS BASED ON DOCUMENTATION IN SUBROUTINE 
C     OUTHYDRO OF THE NGM:
C     
C     THE FUNDAMENTAL HYDROSTATIC EQUATION IS
C        D(HEIGHT)
C        ---------  =  TAU = VIRTUAL TEMPERATURE * (RGAS/GRAVITY)
C        D (Z)
C      WHERE
C        Z = MINUS LOG OF PRESSURE (-LN(P)).
C
C     SEA-LEVEL PRESSURE IS COMPUTED FROM THE FORMULA
C        PRESS(MSL) = PRESS(GROUND) * EXP( F)
C     WHERE
C        F        = HEIGHT OF GROUND / MEAN TAU
C        MEAN TAU = ( TAU(GRND) + TAU(SL) ) / 2
C     
C     IN THE NGM TAU(GRND) AND TAU(SL) ARE FIRST SET USING A 
C     6.5DEG/KM LAPSE RATE FROM SIGMA LEVEL K=1.  THIS IS MODIFIED
C     BY A CORRECTION BASED ON THE CRITICAL TAU OF THE SHUELL
C     CORRECTION:
C                  TAUCR=(RGASD/GRAVITY) * 290.66
C   
C     1) WHERE ONLY TAU(SL) EXCEEDS TAUCR, CHANGE TAU(SL) TO TAUCR.
C
C     2) WHERE BOTH TAU(SL) AND TAU(GRND) EXCEED TAUCR,
C        CHANGE TAU(SL) TO TAUCR-CONST*(TAU(GRND)-TAUCR  )**2
C        WHERE CONST = .005 (GRAVITY/RGASD)
C   
C     THE AVERAGE OF TAU(SL) AND TAU(GRND) IS THEN USED TOGETHER
C     WITH THE GROUND HEIGHT AND PRESSURE TO DERIVE THE PRESSURE
C     AT SEA LEVEL. 
C     
C     EXTRAPOLATING FROM TEMPERATURES IN THE FIRST ATMOSPHERIC ETA
C     LAYER LEAD TO EXCESSIVELY NOISY UNDERGROUND TEMPERATURE FIELDS.
C     AFTER EXPERIMENTATION WE OPTED TO USE A MEAN TEMPERATURE FROM
C     THE LOWEST DP=30MB ABOVE THE GROUND.  FROM THIS MEAN TEMPERATURE
C     WE EXTRAPOLATE TO GET THE GROUND AND SEA LEVEL TEMPERATURES.
C     USING A TRUE MASS WEIGHTED LAYER MEAN INSTEAD OF THE SIMPLE
C     ARITHMETIC MEAN CODED SHOWED LITTLE QUALITATIVE DIFFERENCE.
C     
C     HEIGHT OF THE 1000MB SURFACE IS COMPUTED FROM THE MSL PRESSURE
C     FIELD USING THE FORMULA:
C     
C       P(MSL) - P(1000MB) = MEAN DENSITY * GRAVITY * HGT(1000MBS)
C     
C     WHERE P(MSL) IS THE SEA LEVEL PRESSURE FIELD WE HAVE JUST
C     COMPUTED.
C     
C     BELOW GROUND VIRTUAL TEMPERATURES ARE COMPUTED ASSUMING A CONSTANT
C     6.5DEG/KM LAPSE RATE.  TO CONVERT VIRTUAL TEMPERATURE TO A DRY
C     TEMPERATURE WE USE THE DP=30MB MEAN SPECIFIC HUMIDITY.
C     
C     LINE COMMENTED OUT WITH 'CX' ARE LEFT OVER FROM EARLIER
C     EXPERIMENTATION.  THEY REMAIN TO REMIND READERS WHAT 
C     WAS TRIED.
C   .     
C     
C PROGRAM HISTORY LOG:
C   93-02-02  RUSS TREADON
C   98-06-08  T BLACK - CONVERSION FROM 1-D TO 2-D
C   00-01-04  JIM TUCCILLO - MPI VERSION
C     
C USAGE:    CALL NGMSLP2
C   INPUT ARGUMENT LIST:
C     NONE     
C
C   OUTPUT ARGUMENT LIST: 
C     NONE
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UTILITIES:
C       NONE
C     LIBRARY:
C       COMMON   - VRBLS
C                  EXTRA
C                  MAPOT
C                  LOOPS
C                  DYNAMD
C                  MASKS
C                  PVRBLS
C                  INDX
C     
C   ATTRIBUTES:
C     LANGUAGE: FORTRAN
C     MACHINE : CRAY C-90
C$$$  
C     
C     
C     INCLUDE PARAMETERS.  SET LOCAL PARAMETERS.
      INCLUDE "parmeta"
      INCLUDE "params"
      PARAMETER (GAMMA=6.5/1000.,ZSL=0.0)
      PARAMETER (TAUCR=RD*GI*290.66,CONST=0.005*G/RD)
      PARAMETER (GORD=G/RD,DP=60.E2)
      PARAMETER (ISMTHT=4,ISMTHQ=2)
C     
C     DECLARE VARIABLES
      REAL TBAR(IM,JM),QBAR(IM,JM),ALPBAR(IM,JM)
      REAL TSFC(IM,JM),QSFC(IM,JM),EGRID(IM,JM)
C     
C     INCLUDE COMMON BLOCKS
      INCLUDE "VRBLS.comm"
      INCLUDE "EXTRA.comm"
      INCLUDE "MAPOT.comm"
      INCLUDE "LOOPS.comm"
      INCLUDE "DYNAMD.comm"
      INCLUDE "MASKS.comm"
      INCLUDE "PVRBLS.comm"
      INCLUDE "INDX.comm"
      INCLUDE "CTLBLK.comm"
C     
C**********************************************************************
C     START NGMSLP HERE.
C     
C     LOOP OVER HORIZONTAL GRID.
C
!$omp  parallel do
!$omp& private(alpm,alps,l,llmh,nlev,pm,psfc,ptop,qsum,thsfc,tsum)
      doout30: DO J=JSTA,JEND
      doin30: DO I=1,IM
         LLMH = LMH(I,J)
C     
C        LOCATE TOP OF LAYER OVER WHICH TO COMPUTE MEAN FIELDS.
         PSFC    = PD(I,J)+PT
         PTOP    = PSFC-DP
         QSFC(I,J) = (1.-SM(I,J))*QS(I,J)+SM(I,J)*QZ0(I,J)
         QSFC(I,J) = AMAX1(H1M12,QSFC(I,J))
         THSFC   = (1.-SM(I,J))*THS(I,J)+SM(I,J)*THZ0(I,J)
         TSFC(I,J) = THSFC*(P1000/PSFC)**(-CAPA)
C
C        COMPUTE MEAN FIELDS.
         NLEV = 1
         ALPS = D50*(ALPINT(I,J,LLMH)+ALPINT(I,J,LLMH+1))
         TSUM = T(I,J,LLMH)
         QSUM = Q(I,J,LLMH)
         TBAR(I,J)   = TSUM
         QBAR(I,J)   = QSUM
         ALPBAR(I,J) = ALPS
         IF (LLMH.EQ.LM) CYCLE doin30
         do10: DO L = LLMH-1,1,-1
            ALPM = D50*(ALPINT(I,J,L)+ALPINT(I,J,L+1))
            PM   = EXP(ALPM)
            IF (PM.LT.PTOP) EXIT do10 
            NLEV = NLEV + 1
            ALPS = ALPS + ALPM
            TSUM = TSUM + T(I,J,L)
            QSUM = QSUM + Q(I,J,L)
         END DO do10
         RNLEV     = 1./NLEV
         ALPBAR(I,J) = ALPS*RNLEV
         TBAR(I,J)   = TSUM*RNLEV
         QBAR(I,J)   = QSUM*RNLEV
         END DO doin30
         END DO doout30
C     
C     SMOOTH LAYER MEAN TEMPERATURE AND SPECIFIC HUMIDITY.
C     BOUND SPECIFIC HUMIDITY SO THAT IT IS NON-NEGATIVE.
C     
      CALL P2FILT(ISMTHT,HBM2,TBAR)
      CALL P2FILT(ISMTHT,HBM2,TSFC)
      CALL P2FILT(ISMTHQ,HBM2,QBAR)
      CALL P2FILT(ISMTHQ,HBM2,QSFC)
      CALL BOUNDL(QBAR,H1M12,H99999,IM,JM)
      CALL BOUNDL(QSFC,H1M12,H99999,IM,JM)
C     
C     LOOP OVER HORIZONTAL GRID.
C     
!$omp  parallel do
!$omp& private(alpavg,alpsfc,llmh,pavg,psfc,ptop,qavg,rhoavg,rrhog,
!$omp&         tau,tauavg,tausfc,tausl,tavg,tvrbar,tvrsfc,tvrsl,
!$omp&         tvrt,tvrtal,zbar,zl,zsfc)
       doout50: DO J=JSTA,JEND
       doin50: DO I=1,IM
         LLMH = LMH(I,J)
         ZSFC = FIS(I,J)*GI
         PSFC = PD(I,J)+PT
         ALPSFC = ALOG(PSFC)
         PTOP   = PSFC-DP
         SLP(I,J) = PSFC
CX         IF (LLMH.EQ.LM) GOTO 50
C
C        COMPUTE LAYER MEAN TAU (VIRTUAL TEMP*RD/G).
         TVRT = TBAR(I,J)*(H1+D608*QBAR(I,J))
         TAU  = TVRT*RD*GI
C
C        COMPUTE HEIGHT OF MEAN ATMOSPHERIC LAYER.
CX         QSFC   = (1.-SM(I,J))*QS(I,J)+SM(I,J)*QZ0(I,J)
CX         QSFC   = AMAX1(H1M12,QSFC)
         QAVG   = D50*(QSFC(I,J)+QBAR(I,J))
CX         THSFC  = (1.-SM(I,J))*THS(I,J)+SM(I,J)*THZ0(I,J)
CX         TSFC   = THSFC*(P1000/PSFC)**(-CAPA)
         TAVG   = D50*(TSFC(I,J)+TBAR(I,J))
         TVRBAR = TAVG*(H1+D608*QAVG)
         ZBAR   = RD*GI*TVRBAR*(ALPSFC-ALPBAR(I,J)) + ZSFC
C     
C        COMPUTE TAU AT THE GROUND (Z=ZSFC) AND SEA LEVEL (Z=0)
C        ASSUMING A CONSTANT LAPSE RATE OF GAMMA=6.5DEG/KM.
         TVRSFC = TVRT + (ZBAR-ZSFC)*GAMMA
         TAUSFC = TVRSFC*RD*GI
         TVRSL  = TVRT + (ZBAR- ZSL)*GAMMA
         TAUSL  = TVRSL*RD*GI
C     
C        IF NEED BE APPLY SHEULL CORRECTION.
         IF ((TAUSL.GT.TAUCR).AND.(TAUSFC.LE.TAUCR)) THEN
            TAUSL=TAUCR
         ELSEIF ((TAUSL.GT.TAUCR).AND.(TAUSFC.GT.TAUCR)) THEN
            TAUSL = TAUCR-CONST*(TAUSFC-TAUCR)**2
         ENDIF
C     
C        COMPUTE MEAN TAU.
         TAUAVG = D50*(TAUSL+TAUSFC)
C     
C        COMPUTE SEA LEVEL PRESSURE.
         IF (LLMH.LT.LM) SLP(I,J) = PSFC*EXP(ZSFC/TAUAVG)
C     
C        COMPUTE 1000MB HEIGHTS.
         ALPAVG   = D50*(ALOG(PSFC)+ALOG(SLP(I,J)))
         PAVG     = EXP(ALPAVG)
         RHOAVG   = PAVG*GI/TAUAVG
         RRHOG    = H1/(RHOAVG*G)
         Z1000(I,J) = (SLP(I,J)-P1000)*RRHOG
C     
C        COMPUTE TEMPERATURES ON ETA LEVELS BELOW GROUND.
         IF ((ZSFC-ZSL).EQ.0.) CYCLE doin50
CX       TLAPSE = GORD*(TAUSFC-TAUSL)/(ZSFC-ZSL)
         IF (LLMH.EQ.LM) CYCLE doin50
         do40: DO L = LLMH+1,LM
            ZL     = GI*DFL(L)
            TVRTL  = TVRT + (ZBAR-ZL)*GAMMA
            T(I,J,L) = TVRTL/(H1+D608*QBAR(I,J))
         END DO do40
C     
C     MOVE TO NEXT HORIZONTAL GRIDPOINT.
      END DO doin50
      END DO doout50
C     
C     WE NOW HAVE THE SHEULL SEA LEVEL PRESSURE FIELD
C     IN ARRAY SLP.  APPLY A FIVE POINT SMOOTHER.
C     
!$omp  parallel do
      DO J=JSTA,JEND
      DO I=1,IM
        EGRID(I,J)=SLP(I,J)
      ENDDO
      ENDDO
C
      CALL EXCH(SLP) 
      do90: DO KS=1,KSLPD
         do80: DO J=JSTA_M2,JEND_M2
         IEND=IM-1-MOD(J+1,2)
         DO I=2,IEND
           IF (FIS(I,J).GT.H1000)         GO TO 70
           IF (FIS(I+IHW(J),J-1).GT.H1000)GO TO 70
           IF (FIS(I+IHE(J),J-1).GT.H1000)GO TO 70
           IF (FIS(I+IHW(J),J+1).GT.H1000)GO TO 70
           IF (FIS(I+IHE(J),J+1).GT.H1000)GO TO 70
           IF (KS.GT.1) CYCLE do80
 70        EGRID(I,J)=D125*
     1          (SLP(I+IHW(J),J-1)+SLP(I+IHE(J),J-1)
     1          +SLP(I+IHW(J),J+1)+SLP(I+IHE(J),J+1)
     2          +H4*SLP(I,J))
         ENDDO
         END DO do80
         DO J=JSTA,JEND
         DO I=1,IM
           SLP(I,J) = EGRID(I,J)
         ENDDO
         ENDDO
      END DO do90
C     
C     END OF ROUTINE.
C     

      RETURN
      END
