MAPSSLP.f

!
      SUBROUTINE mapsslp(TPRES)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
!   INPUT ARGUMENT LIST:
!     TPRES    - TEMPERATURE at pressure levels 
!
!   OUTPUT ARGUMENT LIST:
!     PSLP - THE FINAL REDUCED SEA LEVEL PRESSURE ARRAY
!
!-----------------------------------------------------------------------
      use ctlblk_mod, only: jsta, jend, spl, smflag, lm, im, jsta_2l, jend_2u, &
                            lsm, jm, grib, spval
      use gridspec_mod, only: maptype, dxval
      use vrbls3d, only: pmid, t, pint
      use vrbls2d, only: pslp, fis
      use masks, only: lmh
      use params_mod, only: rog, p1000, capa, erad, pi ,gi

      implicit none
!      
      include "mpif.h"
!
      REAL tpres(im,jsta_2l:jend_2u,lsm)
 
      real lapses, expo,expinv,tsfcnew

      REAL,dimension(im, jsta_2l:jend_2u) :: t700
      real,dimension(im,2) :: sdummy
      REAL,dimension(im,jm) :: grid1, th700
      INTEGER nsmooth
      integer l, j, i, k, ii, jj 
      real dxm
!-----------------------------------------------------------------------
!***
        lapses = 0.0065
! deg K / meter
        expo   = rog*lapses
        expinv = 1./expo

      DO l=1,lsm

!$omp parallel do private(i,j)
        DO j=jsta,jend
          DO i=1,im
            if(spl(l) == 70000.)THEN
              if(tpres(i,j,l) <spval)THEN
                t700(i,j)  = tpres(i,j,l) 
                th700(i,j) = t700(i,j)*(p1000/70000.)**capa
              else
                t700(i,j)  = spval
                th700(i,j) = spval
              endif
            endif
          ENDDO
        ENDDO

      ENDDO


! smooth 700 mb temperature first
       if(maptype==6) then
         if(grib=='grib2') then
            dxm=(dxval / 360.)*(erad*2.*pi)/1.d6  ! [mm]
         endif
       else
         dxm = dxval
       endif
       if(grib == 'grib2')then
         dxm = dxm/1000.0 ! [m]
       endif

       IF (smflag) THEN   
         nsmooth=nint(10.*(13500./dxm))
        call allgetherv(th700)
        do k = 1,nsmooth
          CALL smooth(th700,sdummy,im,jm,0.5)
        end do
        ENDIF
          ii=im/2
          jj=(jsta+jend)/2
!          if(i==ii.and.j==jj)                              &
!             print*,'Debug TH700(i,j), i,j',TH700(i,j), i,j

       DO j=jsta,jend
         DO i=1,im
         if(t700(i,j) <spval) then
         t700(i,j) = th700(i,j)*(70000./p1000)**capa
          IF (t700(i,j)>100.) THEN
           tsfcnew = t700(i,j)*(pmid(i,j,lm)/70000.)**expo
!     effective sfc T based on 700 mb temp
          ELSE
           tsfcnew = t(i,j,lm)
           ENDIF
          pslp(i,j) = pint(i,j,nint(lmh(i,j))+1)*               &
              ((tsfcnew+lapses*fis(i,j)*gi)/tsfcnew)**expinv
!          print*,'PSLP(I,J),I,J',PSLP(I,J),I,J
           grid1(i,j)=pslp(i,j)
         else
           pslp(i,j) = spval
           grid1(i,j) = spval
         endif

         ENDDO
       ENDDO

         IF (smflag) THEN
! - in WRF number of passes depends on the resolution: nsmooth=int(15*(13/dxval))
         nsmooth=nint(15.*(13500./dxm))
         call allgetherv(grid1)
         do k=1,nsmooth
          CALL smooth(grid1,sdummy,im,jm,0.5)
         end do
         DO j=jsta,jend
            DO i=1,im
               pslp(i,j)=grid1(i,j)
            ENDDO
         ENDDO
         ENDIF
!

      RETURN
      END