SLP_NMM.f

      SUBROUTINE memslp_nmm(TPRES,QPRES,FIPRES)
!
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
!   SUBROUTINE:  MEMSLP      MEMBRANE SLP REDUCTION
!
! ABSTRACT:  THIS ROUTINE COMPUTES THE SEA LEVEL PRESSURE
!            REDUCTION USING THE MESINGER RELAXATION
!            METHOD FOR SIGMA COORDINATES.
!            A BY-PRODUCT IS THE
!            SET OF VALUES FOR THE UNDERGROUND TEMPERATURES
!            ON THE SPECIFIED PRESSURE LEVELS
!
! PROGRAM HISTORY LOG:
!   99-09-23  T BLACK - REWRITTEN FROM ROUTINE SLP (ETA
!                       COORDINATES)
!   02-07-26  H CHUANG - PARALLIZE AND MODIFIED FOR WRF A/C GRIDS
!                        ALSO REDUCE S.O.R. COEFF FROM 1.75 to 1.25 
!                        BECAUSE THERE WAS NUMERICAL INSTABILITY  
!   02-08-21  H CHUANG - MODIFIED TO ALWAYS USE OLD TTV FOR RELAXATION
!                        SO THAT THERE WAS BIT REPRODUCIBILITY BETWEEN
!                        USING ONE AND MULTIPLE TASKS         
!   13-12-06  H CHUANG - REMOVE EXTRA SMOOTHING OF SLP AT THE END
!   19-10-30  Bo CUI - REMOVE "GOTO" STATEMENT
!
! USAGE:  CALL SLPSIG FROM SUBROUITNE ETA2P
!
!   INPUT ARGUMENT LIST:
!     PD   - SFC PRESSURE MINUS PTOP
!     FIS  - SURFACE GEOPOTENTIAL
!     T    - TEMPERATURE 
!     Q    - SPECIFIC HUMIDITY
!     FI   - GEOPOTENTIAL
!     PT   - TOP PRESSURE OF DOMAIN
!
!   OUTPUT ARGUMENT LIST:
!     PSLP - THE FINAL REDUCED SEA LEVEL PRESSURE ARRAY
!
!   SUBPROGRAMS CALLED:
!     UNIQUE:
!             NONE
!
!-----------------------------------------------------------------------
      use vrbls3d, only: pint, zint, t, q
      use vrbls2d, only: pslp, fis
      use masks, only: lmh
      use params_mod, only: overrc, ad05, cft0, g, rd, d608, h1, kslpd
      use ctlblk_mod, only: jsta, jend, spl, num_procs, mpi_comm_comp, lsmp1, jsta_m2, jend_m2,&
              lm, jsta_m, jend_m, im, jsta_2l, jend_2u, im_jm, lsm, jm
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
      include "mpif.h"
!-----------------------------------------------------------------------
      integer, PARAMETER :: nfill=0,nrlx1=500,nrlx2=100
!-----------------------------------------------------------------------
      real,dimension(IM,JSTA_2L:JEND_2U,LSM),intent(in) :: qpres
      real,dimension(IM,JSTA_2L:JEND_2U,LSM),intent(inout) :: tpres,fipres
      REAL :: ttv(im,jsta_2l:jend_2u),tnew(im,jsta_2l:jend_2u)    & 
              ,slpx(im,jsta_2l:jend_2u)                           &
              ,p1(im,jsta_2l:jend_2u),htm2d(im,jsta_2l:jend_2u)
      REAL :: htmo(im,jsta_2l:jend_2u,lsm)              
      real p2,gz1,gz2,tlyr,spll,pchk,psfc,slope,tvrt,dis,tinit
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
      INTEGER :: kmntm(lsm),imnt(im_jm,lsm),jmnt(im_jm,lsm)        &
                 ,lmho(im,jsta_2l:jend_2u)
      INTEGER :: ihe(jm),ihw(jm),ive(jm),ivw(jm),ihs(jm),ihn(jm)
      integer ii,jj,i,j,l,n,km,ks,kp,kmn,kmm,kount,lp,llmh,lhmnt   &
              ,lmhij,lmap1,lxxx,ierr,nrlx,ihh2
!-----------------------------------------------------------------------
      LOGICAL :: done(im,jsta_2l:jend_2u)
      logical, parameter :: debugprint = .false.
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!***
!***  CALCULATE THE I-INDEX EAST-WEST INCREMENTS
!***
!
      ii=279
      jj=314
      DO j=1,jm
        ihe(j)=mod(j+1,2)
        ihw(j)=ihe(j)-1
      ENDDO
!      print*,'relaxation coeff= ',OVERRC
!-----------------------------------------------------------------------
!***
!***  INITIALIZE ARRAYS.  LOAD SLP ARRAY WITH SURFACE PRESSURE.
!***
!$omp parallel do 
      DO j=jsta,jend
      DO i=1,im
        llmh=nint(lmh(i,j))
        pslp(i,j)=pint(i,j,llmh+1)
        if(debugprint .and. i==ii .and. j==jj)print*,'Debug: FIS,IC for PSLP='      &
        ,fis(i,j),pslp(i,j)
        ttv(i,j)=0.
        lmho(i,j)=0
    done(i,j)=.false.
      ENDDO
      ENDDO
!
!***  CALCULATE SEA LEVEL PRESSURE FOR PROFILES (AND POSSIBLY
!***  FOR POSTING BY POST PROCESSOR).
!
!--------------------------------------------------------------------
!***
!***  CREATE A 3-D "HEIGHT MASK" FOR THE SPECIFIED PRESSURE LEVELS
!***  (1 => ABOVE GROUND) AND A 2-D INDICATOR ARRAY THAT SAYS 
!***  WHICH PRESSURE LEVEL IS THE LOWEST ONE ABOVE THE GROUND
!***
      DO 100 l=1,lsm
      spll=spl(l)
!       
      DO j=jsta,jend
      DO i=1,im
        psfc=pslp(i,j)
        pchk=psfc
        IF(nfill>0)THEN
     pchk=pint(i,j,nint(lmh(i,j))+1-nfill)
        ENDIF
!        IF(SM(I,J)>0.5.AND.FIS(I,J)<1.)PCHK=PSLP(I,J)
        IF(fis(i,j)<1.)pchk=pslp(i,j)
!
!       IF(SPLL<PSFC)THEN
        IF(spll<pchk)THEN
          htmo(i,j,l)=1.
        ELSE
          htmo(i,j,l)=0.
          IF(l>1.AND.htmo(i,j,l-1)>0.5)lmho(i,j)=l-1
        ENDIF
!
        IF(l==lsm.AND.htmo(i,j,l)>0.5)lmho(i,j)=lsm
        if(debugprint .and. i==ii .and. j==jj)print*,'Debug: HTMO= ',htmo(i,j,l)
      ENDDO
      ENDDO
!
  100 CONTINUE
!      if(jj>=jsta.and.jj<=jend)
!     +print*,'Debug: LMHO=',LMHO(ii,jj)
!--------------------------------------------------------------------
!***
!***  WE REACH THIS LINE IF WE WANT THE MESINGER ETA SLP REDUCTION
!***  BASED ON RELAXATION TEMPERATURES.  THE FIRST STEP IS TO
!***  FIND THE HIGHEST LAYER CONTAINING MOUNTAINS.
!***
      loop210: DO l=lsm,1,-1
!
      DO j=jsta,jend
      DO i=1,im
        IF(htmo(i,j,l)<0.5) cycle loop210
      ENDDO
      ENDDO
!
      lhmnt=l+1
      exit loop210
      enddo loop210

      if(debugprint)print*,'Debug in SLP: LHMNT=',lhmnt
      if ( num_procs > 1 ) then
      CALL mpi_allreduce                                          &  
       (lhmnt,lxxx,1,mpi_integer,mpi_min,mpi_comm_comp,ierr)
      lhmnt = lxxx
      end if
    
      IF(lhmnt==lsmp1)THEN
        go to 325
      ENDIF
      if(debugprint)print*,'Debug in SLP: LHMNT A ALLREDUCE=',lhmnt
!***
!***  NOW GATHER THE ADDRESSES OF ALL THE UNDERGROUND POINTS.
!***
!$omp parallel do private(kmn,kount)
      DO 250 l=lhmnt,lsm
      kmn=0
      kmntm(l)=0
      kount=0
      DO 240 j=jsta_m2,jend_m2
!      DO 240 J=JSTA_M,JEND_M
      DO 240 i=2,im-1
      kount=kount+1
      imnt(kount,l)=0
      jmnt(kount,l)=0
      IF(htmo(i,j,l)>0.5) cycle             
      kmn=kmn+1
      imnt(kmn,l)=i
      jmnt(kmn,l)=j
  240 CONTINUE
      kmntm(l)=kmn
  250 CONTINUE
!
!
!***  CREATE A TEMPORARY TV ARRAY, AND FOLLOW BY SEQUENTIAL
!***  OVERRELAXATION, DOING NRLX PASSES.
!
!     IF(NTSD==1)THEN
        nrlx=nrlx1
!     ELSE
!       NRLX=NRLX2
!     ENDIF
!
!!$omp parallel do private(i,j,tinit,ttv)
      DO 300 l=lhmnt,lsm
!
      DO 270 j=jsta,jend
      DO 270 i=1,im
      ttv(i,j)=tpres(i,j,l)
      IF(ttv(i,j)<150. .and. ttv(i,j)>325.0)print*            &  
        ,'abnormal IC for T relaxation',i,j,ttv(i,j)
      htm2d(i,j)=htmo(i,j,l)
  270 CONTINUE
!
!***  FOR GRID BOXES NEXT TO MOUNTAINS, COMPUTE TV TO USE AS
!***  BOUNDARY CONDITIONS FOR THE RELAXATION UNDERGROUND
!
      CALL exch2(htm2d(1,jsta_2l))   !NEED TO EXCHANGE TWO ROW FOR E GRID
      DO j=jsta_m2,jend_m2
      DO i=2,im-1
        IF(htm2d(i,j)>0.5.AND.htm2d(i+ihw(j),j-1)*htm2d(i+ihe(j),j-1) &
          *htm2d(i+ihw(j),j+1)*htm2d(i+ihe(j),j+1)                       &
          *htm2d(i-1     ,j  )*htm2d(i+1     ,j  )                       &
          *htm2d(i       ,j-2)*htm2d(i       ,j+2)<0.5)THEN
!HC MODIFICATION FOR C AND A GRIDS
!HC        IF(HTM2D(I,J)>0.5.AND.
!HC     1     HTM2D(I-1,J)*HTM2D(I+1,J)
!HC     2    *HTM2D(I,J-1)*HTM2D(I,J+1)
!HC     3    *HTM2D(I-1,J-1)*HTM2D(I+1,J-1)
!HC     4    *HTM2D(I-1,J+1)*HTM2D(I+1,J+1)<0.5)THEN
!     
          ttv(i,j)=tpres(i,j,l)*(1.+0.608*qpres(i,j,l))
        ENDIF
!       if(i==ii.and.j==jj)print*,'Debug:L,TTV B SMOO= ',l,TTV(I,J) 
      ENDDO
      ENDDO
!
      kmm=kmntm(l)
!
      DO 285 n=1,nrlx
      CALL exch2(ttv(1,jsta_2l))
!      print*,'Debug:L,KMM=',L,KMM
      DO 280 km=1,kmm
      i=imnt(km,l)
      j=jmnt(km,l)
      tinit=ttv(i,j)
      tnew(i,j)=ad05*(4.*(ttv(i+ihw(j),j-1)+ttv(i+ihe(j),j-1)            &
                        +ttv(i+ihw(j),j+1)+ttv(i+ihe(j),j+1))            &
                        +ttv(i-1,j)       +ttv(i+1,j)                    &
                        +ttv(i,j-2)       +ttv(i,j+2))                   &
                        -cft0*ttv(i,j)
!HC MODIFICATION FOR C AND A GRIDS
! eight point relaxation using old TTV
!HC       TNEW(I,J)=AD05*(4.*(TTV(I-1,J)+TTV(I+1,J)
!HC     1                  +TTV(I,J-1)+TTV(I,J+1))
!HC     2                  +TTV(I-1,J-1)+TTV(I+1,J-1)
!HC     3                  +TTV(I-1,J+1)+TTV(I+1,J+1))
!HC     4                  -CFT0*TTV(I,J)
!
!     if(i==ii.and.j==jj)print*,'Debug: L,TTV A S'
!    1,l,TTV(I,J),N
!     1,l,TNEW(I,J),N
  280 CONTINUE
!
      DO km=1,kmm
       i=imnt(km,l)
       j=jmnt(km,l)
       ttv(i,j)=tnew(i,j)
      END DO
  285 CONTINUE
!
      DO 290 km=1,kmm
      i=imnt(km,l)
      j=jmnt(km,l)
      tpres(i,j,l)=ttv(i,j)
  290 CONTINUE
  300 CONTINUE
!----------------------------------------------------------------
!***
!***  CALCULATE THE SEA LEVEL PRESSURE AS PER THE NEW SCHEME.
!***  INTEGRATE THE HYDROSTATIC EQUATION DOWNWARD FROM THE
!***  GROUND THROUGH EACH OUTPUT PRESSURE LEVEL (WHERE TV
!***  IS NOW KNOWN) TO FIND GZ AT THE NEXT MIDPOINT BETWEEN
!***  PRESSURE LEVELS.  WHEN GZ=0 IS REACHED, SOLVE FOR THE
!***  PRESSURE.
!***
!
!***  COUNT THE POINTS WHERE SLP IS DONE BELOW EACH OUTPUT LEVEL
!
      kount=0
      DO j=jsta,jend
      DO i=1,im
!        P1(I,J)=SPL(NINT(LMH(I,J)))
!        DONE(I,J)=.FALSE.
        IF(abs(fis(i,j))<1.)THEN
          pslp(i,j)=pint(i,j,nint(lmh(i,j))+1)
          done(i,j)=.true.
          kount=kount+1
          if(i==ii.and.j==jj)print*,'Debug:DONE,PSLP A S1='       &  
               ,done(i,j),pslp(i,j)
        ELSE IF(fis(i,j)<-1.0) THEN
          DO l=lm,1,-1
            IF(zint(i,j,l)>0.)THEN
              pslp(i,j)=pint(i,j,l)/exp(-zint(i,j,l)*g                &
              /(rd*t(i,j,l)*(q(i,j,l)*d608+1.0)))
              done(i,j)=.true.
              if(debugprint .and. i==ii.and.j==jj)print*                           &
              ,'Debug:DONE,PINT,PSLP A S1='                           &
               ,done(i,j),pint(i,j,l),pslp(i,j)
              EXIT 
            END IF
          END DO
        ENDIF
      ENDDO
      ENDDO
!
      kmm=kmntm(lsm)
!$omp parallel do private(gz1,gz2,i,j,lmap1,p1,p2),shared(pslp)

loop320: DO km=1,kmm
      i=imnt(km,lsm)
      j=jmnt(km,lsm)
      IF(done(i,j)) cycle 
      lmhij=lmho(i,j)
      gz1=fipres(i,j,lmhij)
      p1(i,j)=spl(lmhij)
!
      lmap1=lmhij+1
      DO l=lmap1,lsm
        p2=spl(l)
        tlyr=0.5*(tpres(i,j,l)+tpres(i,j,l-1))
        gz2=gz1+rd*tlyr*alog(p1(i,j)/p2)
        fipres(i,j,l)=gz2
!        if(i==ii.and.j==jj)print*,'Debug:L,FI A S2=',L,GZ2
        IF(gz2<=0.)THEN
          pslp(i,j)=p1(i,j)/exp(-gz1/(rd*tpres(i,j,l-1)))
!          if(i==ii.and.j==jj)print*,'Debug:PSLP A S2=',PSLP(I,J)
          done(i,j)=.true.
          kount=kount+1
          cycle loop320
        ENDIF
        p1(i,j)=p2
        gz1=gz2
      ENDDO
!HC EXPERIMENT
      lp=lsm
      slope=-6.6e-4 
      tlyr=tpres(i,j,lp)-0.5*fipres(i,j,lp)*slope
      pslp(i,j)=spl(lp)/exp(-fipres(i,j,lp)/(rd*tlyr))
      done(i,j)=.true.
!      if(i==ii.and.j==jj)print*,'Debug:spl,FI,TLYR,PSLPA3='    &  
!         ,spl(lp),FIPRES(I,J,LP),TLYR,PSLP(I,J)       
!HC EXPERIMENT
ENDDO loop320
!
!***  WHEN SEA LEVEL IS BELOW THE LOWEST OUTPUT PRESSURE LEVEL,
!***  SOLVE THE HYDROSTATIC EQUATION BY CHOOSING A TEMPERATURE
!***  AT THE MIDPOINT OF THE LAYER BETWEEN THAT LOWEST PRESSURE
!***  LEVEL AND THE GROUND BY EXTRAPOLATING DOWNWARD FROM T ON
!***  THE LOWEST PRESSURE LEVEL USING THE DT/DFI BETWEEN THE
!***  LOWEST PRESSURE LEVEL AND THE ONE ABOVE IT.
!
!      TOTAL=(IM-2)*(JM-4)
!
!HC      DO 340 LP=LSM,1,-1
!      IF(KOUNT==TOTAL)GO TO 350
!HC MODIFICATION FOR SMALL HILL HIGH PRESSURE SITUATION
!HC IF SURFACE PRESSURE IS CLOSER TO SEA LEVEL THAN LWOEST
!HC OUTPUT PRESSURE LEVEL, USE SURFACE PRESSURE TO DO EXTRAPOLATION
 325  CONTINUE
      lp=lsm
      DO 330 j=jsta,jend
      DO 330 i=1,im
      if(debugprint .and. i==ii.and.j==jj)print*,'Debug: with 330 loop'
      IF(done(i,j)) cycle   
      if(debugprint .and. i==ii.and.j==jj)print*,'Debug: still within 330 loop'
!HC Comment out the following line for situation with terrain 
!HC at boundary (ie FIPRES<0)
!HC because they were not counted as undergound point for 8 pt
!HC relaxation
!HC      IF(FIPRES(I,J,LP)<0.)GO TO 330
!      IF(FIPRES(I,J,LP)<0.)THEN  
!       DO LP=LSM,1,-1
!        IF (FIPRES(I,J) <= 0)

!      IF(FIPRES(I,J,LP)<0..OR.DONE(I,J))GO TO 330
!     SLOPE=(TPRES(I,J,LP)-TPRES(I,J,LP-1))
!     & /(FIPRES(I,J,LP)-FIPRES(I,J,LP-1))     
      slope=-6.6e-4
      IF(pint(i,j,nint(lmh(i,j))+1)>spl(lp))THEN
       llmh=nint(lmh(i,j))
       tvrt=t(i,j,llmh)*(h1+d608*q(i,j,llmh))
       dis=zint(i,j,llmh+1)-zint(i,j,llmh)+0.5*zint(i,j,llmh+1)
       tlyr=tvrt-dis*g*slope
       pslp(i,j)=pint(i,j,llmh+1)*exp(zint(i,j,llmh+1)*g/(rd*tlyr))
!       if(i==ii.and.j==jj)print*,'Debug:PSFC,zsfc,TLYR,PSLPA3='
!     1,PINT(I,J,LLMH+1),ZINT(I,J,LLMH+1),TLYR,PSLP(I,J)
      ELSE
       tlyr=tpres(i,j,lp)-0.5*fipres(i,j,lp)*slope
       pslp(i,j)=spl(lp)/exp(-fipres(i,j,lp)/(rd*tlyr))                   
       if(debugprint .and. i==ii.and.j==jj)print*,'Debug:spl,FI,TLYR,PSLPA3='       &
         ,spl(lp),fipres(i,j,lp),tlyr,pslp(i,j)
      END IF
      done(i,j)=.true.
      kount=kount+1
  330 CONTINUE
!HC  340 CONTINUE
!
  350 CONTINUE
!----------------------------------------------------------------
      RETURN
      END