SUBROUTINE BLOSFC
!C$$$  SUBPROGRAM DOCUMENTATION BLOCK
!C                .      .    .     
!C SUBPROGRAM:    BLOSFC2     SETS BELOW SURFACE VALUES
!C   PRGRMMR: TREADON         ORG: W/NP2      DATE: 93-05-07       
!C     
!C ABSTRACT:  THIS ROUTINE SETS BELOW GROUND Q, U, V, 
!C     AND OMEGA.  FOR U, V, AND OMEGA WE SIMPLY FILL 
!C     BELOW GROUND ARRAY ELEMENTS WITH VALUES FROM 
!C     THE FIRST ATMOSPHERIC ETA LAYER (FAL).  FOR Q 
!C     WE FIRST COMPUTE THE FAL RELATIVE HUMIDITY.  USING
!C     THE GIVEN TEMPERATURE AND PRESSURE WE USE THIS 
!C     FAL RH FIELD TO COMPUTE BELOW SURFACE Q WHICH 
!C     MAINTAINS THE FAL RH.
!C   .     
!C     
!C PROGRAM HISTORY LOG:
!C   93-01-27  RUSS TREADON
!C   93-05-07  RUSS TREADON - ADDED DOCBLOC
!C   96-03-07  MIKE BALDWIN - SPEED UP CODE 
!C   98-06-08  T BLACK      - CONVERSION FROM 1-D TO 2-D
!C   98-08-17  MIKE BALDWIN - COMPUTE RH OVER ICE
!C   98-12-22  MIKE BALDWIN - BACK OUT RH OVER ICE
!C   00-01-03  JIM TUCCILLO - MPI VERSION         
!!   06-01-19  Hai Zhang    - Global Eta
!C     
!C USAGE:    CALL BLOSFC2
!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   - MAPOT
!C                  VRBLS
!C                  LOOPS
!C                  EXTRA
!C                  OMGAOT
!C                  MASKS
!C     
!C   ATTRIBUTES:
!C     LANGUAGE: FORTRAN 90
!C     MACHINE : CRAY C-90
!C$$$  
!C     
!C
!C     INCLUDE PARAMETER STATEMENTS.  SET LOCAL PARAMETERS.
       USE mod_vrbls
       USE mod_extra
       USE mod_masks
       USE mod_pvrbls
       include 'param_o.h'
!!JNT       include 'extra.comm'
       include 'mapot.comm'
       include 'const.h'
       include 'dynam.comm'
       include 'params'

      PARAMETER (DPBND=60.E2,ISMTHP=2,ISMTHT=2,ISMTHQ=2,ISMTHR=2,  &
       CLIMIT=1.E-20)
!C     
!C     DECLARE VARIABLES.
      REAL PBND(IM,JM,nm),QBND(IM,JM,nm),RHBND(IM,JM,nm),TBND(IM,JM,nm)
      REAL PSUM(IM,JM,nm),ICEB(IM,JM,nm),IWM1(IM,JM,nm)
!C     
!C********************************************************************
!C     START BLOSFC HERE.
!C     
!C     SET BELOW GROUND OMEGA.

      DO L = 1,LM
         do n=1,nm
         DO J=1,jm
         DO I=1,IM
           LLMH = LMH(I,J,n)
!           IF(L.GT.LLMH) OMGA(I,J,n,L) = OMGA(I,J,n,LLMH)
         ENDDO
         ENDDO
         ENDDO
      ENDDO
!C     
!C     SET BELOW GROUND U AND V WIND COMPONENTS.
!!$omp  parallel do
!!$omp& private(llmv)
      DO L = 1,LM
         do n=1,nm
         DO J=1,jm
         DO I=1,IM
           LLMV = LMV(I,J,n)
           IF (L.GT.LLMV) THEN
             U(I,J,n,L) = U(I,J,n,LLMV)
             V(I,J,n,L) = V(I,J,n,LLMV)
           ENDIF
         ENDDO
         ENDDO
         ENDDO
      ENDDO
!C
!C     LOOP OVER HORIZONTAL.  AT EACH MASS POINT COMPUTE 
!C     LAYER MEAN P, T, AND Q IN A DPBND THICK BOUNDARY
!C     LAYER FROM THE SURFACE UP.
!C
!!$omp  parallel do
      do n=1,nm
      DO J=1,jm
      DO I=1,IM 
        PBND(I,J,n)= PD(I,J,n) + PT - 0.5*DPBND
        PSUM(I,J,n)= 0
        TBND(I,J,n)= 0
        QBND(I,J,n)= 0
        ICEB(I,J,n)= 0
        IWM1(I,J,n)= 0
      ENDDO
      ENDDO
      ENDDO
!$omp  parallel do  private(dp,iwm1,pbot,pm,ptop,riw)
! !$omp  private(dp,iwm1,pbot,pm,ptop,riw)
      DO L = 1,LM
        do n=1,nm
        DO J=1,jm
        DO I=1,IM
          PM = 0.50*(PINT(I,J,n,L)+PINT(I,J,n,L+1))
          PTOP = PBND(I,J,n)-DPBND*0.5
          PBOT = PBND(I,J,n)+DPBND*0.5
!C   COMPUTE IW
          RIW=0.
          IF(W(I,J,n,L).GT.CLIMIT) THEN
             IF(T(I,J,n,L).LT.258.15)THEN
               RIW=1.
             ELSEIF(T(I,J,n,L).GE.273.15)THEN
               RIW=0.
             ELSE
               IF(IWM1(I,J,n).EQ.1.0)RIW=1.
             ENDIF
          ELSE
             RIW=0.
          ENDIF
          IWM1(I,J,n)=RIW
!C   COMPUTE IW
          IF (PM.GT.PTOP.AND.PM.LE.PBOT) THEN
             DP = PINT(I,J,n,L+1)-PINT(I,J,n,L)
             PSUM(I,J,n) = PSUM(I,J,n) + DP
             TBND(I,J,n) = TBND(I,J,n) + T(I,J,n,L)*DP
             QBND(I,J,n) = QBND(I,J,n) + Q(I,J,n,L)*DP
             ICEB(I,J,n) = ICEB(I,J,n) + RIW*DP
          ENDIF
        ENDDO
        ENDDO
        ENDDO
      ENDDO
!C
      do n=1,nm
      DO J=1,jm
      DO I=1,IM
         IF (PSUM(I,J,n).NE.0.) THEN
            RPSUM   = 1./PSUM(I,J,n)
            TBND(I,J,n) = TBND(I,J,n)*RPSUM
            QBND(I,J,n) = QBND(I,J,n)*RPSUM
            ICEB(I,J,n) = ICEB(I,J,n)*RPSUM
            IF (ICEB(I,J,n).LT.0.5) ICEB(I,J,n)=0.
         ELSE
            LLMH=LMH(I,J,n)
            TBND(I,J,n) = T(I,J,n,LLMH)
            QBND(I,J,n) = Q(I,J,n,LLMH)
            ICEB(I,J,n) = IWM1(I,J,n)
         ENDIF
      ENDDO
      ENDDO
      ENDDO
!C     USE BOUNDARY LAYER PRESSURE, TEMPERATURE, AND SPECIFIC
!C     HUMIDITY ARRAYS TO COMPUTE BOUNDARY LAYER RELATIVE
!C     HUMIDITY
!C     
!C     SET BELOW GROUND Q TO PRESERVE BOUNDARY LAYER
!C     RELATIVE HUMIDITY.
!C     
!!$omp  parallel do private(ai,bi,llmh,pm,qi,qint,qs,qw,tm,tmt0,tmt15)
!!$omp  private(ai,bi,llmh,pm,qi,qint,qs,qw,tm,tmt0,tmt15)
      CALL CALRH2(PBND,TBND,QBND,ICEB,RHBND,IM,JM,nm)
      DO L = 1,LM
        do n=1,nm
        DO J=1,jm
        DO I=1,IM
          LLMH=LMH(I,J,n)
          IF(L.GT.LLMH)THEN
            PM=0.50*(PINT(I,J,n,L)+PINT(I,J,n,L+1))
            TM=T(I,J,n,L)
!C     
            TMT0=TM-273.16
            TMT15=AMIN1(TMT0,-15.)
            AI=0.008855
            BI=1.
            IF(TMT0.LT.-20.)THEN
              AI=0.007225
              BI=0.9674
            ENDIF
            QW=PQ0/PM*EXP(A2*(TM-A3)/(TM-A4))
            QI=QW*(BI+AI*AMIN1(TMT0,0.))
            QINT=QW*(1.-0.00032*TMT15*(TMT15+15.))
            IF(TMT0.LT.-15.)THEN
               QS=QI
            ELSEIF(TMT0.GE.0.)THEN
               QS=QINT
            ELSE
               IF(ICEB(I,J,n).GT.0.0) THEN
                 QS=QI
               ELSE
                 QS=QINT
               ENDIF
            ENDIF
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
!C             DELETE THIS LINE TO SWITCH BACK TO RH VS ICE
            QS=QW
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
!C
!C
            Q(I,J,n,L)=RHBND(I,J,n)*QS
	    if(i.eq.30.and.j.eq.101.and.n.eq.3.and.l.eq.36)then
	       print *,q(i,j,n,l),rhbnd(i,j,n),qs
            endif
            Q(I,J,n,L)=AMAX1(1.e-12,Q(I,J,n,L))
          ENDIF
        ENDDO
        ENDDO
        ENDDO
      ENDDO
!C     
!C     END OF ROUTINE
!C     
      RETURN
      END