SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,r_radius,thunderRadius,&
                        t_bk,p_bk,q,h_bk,zh,  &
                        mxst_p,NVARCLD_P,numsao,OI,OJ,OCLD,OWX,Oelvtn,Odist,&
                        cld_cover_3d,cld_type_3d,wthr_type,pcp_type_3d,     &
                        watericemax, kwatericemax)
!
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:  cloudCover_Surface  cloud cover analysis using surface observation
!
!   PRGMMR: Ming Hu          ORG: GSD/AMB        DATE: 2006-10-30
!
! ABSTRACT: 
!  This subroutine determines 3D cloud fractional cover using surface observations
!    Code based on RUC assimilation code (hybfront/hybcloud.f)
!
! PROGRAM HISTORY LOG:
!    2009-01-20  Hu  Add NCO document block
!
!
!   input argument list:
!     mype        - processor ID
!     nlon        - no. of lons on subdomain (buffer points on ends)
!     nlat        - no. of lats on subdomain (buffer points on ends)
!     nsig        - no. of levels
!     r_radius    - influence radius of the cloud observation
!     thunderRadius -
!
!     t_bk        - 3D background potentional temperature (K)
!     p_bk        - 3D background pressure  (hPa)
!     q           - 3D moisture (water vapor mixing ratio)
!     h_bk        - 3D background height  (m)
!     zh          - terrain (m)
!
!     mxst_p      -  maximum observation number
!     NVARCLD_P   -  first dimension of OLCD
!     numsao      -  observation number
!     OI          -  observation x location
!     OJ          -  observation y location
!     OLCD        -  cloud amount, cloud height, visibility
!     OWX         -  weather observation
!     Oelvtn      -  observation elevation
!     Odist       -  distance from the nearest station
!
!   output argument list:
!     cld_cover_3d- 3D cloud cover
!     cld_type_3d - 3D cloud type
!     wthr_type   - 3D weather type
!     pcp_type_3d - 3D weather precipitation type
!
! USAGE:
!   INPUT FILES: 
!
!   OUTPUT FILES:
!
!
! REMARKS:
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90 
!   MACHINE:  Linux cluster (WJET)
!
!$$$
!
!_____________________________________________________________________
!

  use kinds, only: r_single,i_kind,r_kind

  implicit none

  integer(i_kind),intent(in) :: mype
  REAL(r_single), intent(in) :: r_radius
  integer(i_kind),intent(in) :: nlat,nlon,nsig
  real(r_single), intent(in) :: thunderRadius
!
!  surface observation
!
  INTEGER(i_kind),intent(in) :: mxst_p,NVARCLD_P

!  PARAMETER (LSTAID_P=9)

  INTEGER(i_kind),intent(in) :: numsao
  real(r_single), intent(in) :: OI(mxst_p)  ! x location
  real(r_single), intent(in) :: OJ(mxst_p)  ! y location
  INTEGER(i_kind),intent(in) :: OCLD(NVARCLD_P,mxst_p)  ! cloud amount, cloud height,
                                               ! visibility
  CHARACTER*10,   intent(in) :: OWX(mxst_p)    ! weather
  real(r_single), intent(in) :: Oelvtn(mxst_p) ! elevation
  real(r_single), intent(in) :: Odist(mxst_p)  ! distance from the nearest station

!
!  background
!
  real(r_single),intent(in) :: t_bk(nlon,nlat,nsig)  ! temperature
  real(r_single),intent(in) :: p_bk(nlon,nlat,nsig)  ! pressure
  real(r_single),intent(in) :: zh(nlon,nlat)         ! terrain
  real(r_single),intent(in) :: q(nlon,nlat,nsig)     ! moisture, water vapor mixing ratio (kg/kg)
  real(r_single),intent(in) :: h_bk(nlon,nlat,nsig)  ! height
!
  REAL(r_single),intent(in)   :: watericemax(mxst_p)  ! max of background total liquid water in station
  INTEGER(i_kind),intent(in):: kwatericemax(nlon,nlat)  ! lowest level of background total liquid water in grid
!
!  Variables for cloud analysis
!
  real (r_single),intent(inout) :: cld_cover_3d(nlon,nlat,nsig)
  integer(i_kind),intent(inout) :: cld_type_3d(nlon,nlat,nsig)
  integer(i_kind),intent(inout) :: wthr_type(nlon,nlat)
  integer(i_kind),intent(inout) :: pcp_type_3d(nlon,nlat,nsig)
!
!  local
!
  real (r_single) :: vis2qc(nlon,nlat)
  real (r_single) :: cloud_zthick_p
  data  cloud_zthick_p    /300._r_kind/
!
  REAL (r_kind)   :: spval_p
  PARAMETER ( spval_p    =  99999.0_r_kind )

  INTEGER(i_kind) :: i,j,k,k1
  INTEGER(i_kind) :: i1,j1,ic
  INTEGER(i_kind) :: nx_p, ny_p, nztn_p
  INTEGER(i_kind) :: ista
  INTEGER(i_kind) :: ich, iob,job 
  
  REAL(r_kind) :: min_dist, dist
  REAL(r_kind) :: zdiff
  REAL(r_kind) :: zlev_clr,cloud_dz,cl_base_ista,betav
!
!
!
  real(r_single):: tbk_k(nlon,nlat,nsig)
  real(r_single):: cv_bk(nlon,nlat,nsig)
  real(r_single):: z_lcl(nlon,nlat)
  REAL(r_kind)  :: cf_model_base,t_model_base, ht_base
  REAL(r_kind)  :: t_dry_adiabat,t_inversion_strength

  LOGICAL :: l_cf,l_inversion
  LOGICAL :: if_cloud_exist

  integer(i_kind) :: firstcloud,cl_base_broken_k
  real(r_single)  ::    underlim
  integer(i_kind) :: npts_near_clr


!====================================================================
!  Begin
!
!  set constant names consistent with original RUC code
!
   nx_p=nlon
   ny_p=nlat
   nztn_p=nsig

   vis2qc=-9999.0_r_kind
   npts_near_clr=0
!
!
!*****************************************************************
!  analysis of surface/METAR cloud observations
! *****************************************************************

   DO ista=1,numsao
     i1 = int(oi(ista)+0.0001_r_kind) 
     j1 = int(oj(ista)+0.0001_r_kind)
     min_dist =  Odist(ista)

!mh - grid point has the closest cloud station

! -- find out if any precip is present
     do ich=1,1
       if ( owx(ista)(ich:ich+1)=='SH'  ) wthr_type(i1,j1)=16
       if ( owx(ista)(ich:ich+1)=='TH' .and. &
                min_dist < thunderRadius) wthr_type(i1,j1)=1
       if ( owx(ista)(ich:ich+1)=='RA'  ) wthr_type(i1,j1)=11
       if ( owx(ista)(ich:ich+1)=='SN'  ) wthr_type(i1,j1)=12
       if ( owx(ista)(ich:ich+1)=='PL'  ) wthr_type(i1,j1)=13
       if ( owx(ista)(ich:ich+1)=='DZ'  ) wthr_type(i1,j1)=14
       if ( owx(ista)(ich:ich+1)=='UP'  ) wthr_type(i1,j1)=15
       if ( owx(ista)(ich:ich+1)=='BR'  ) wthr_type(i1,j1)=21
       if ( owx(ista)(ich:ich+1)=='FG'  ) wthr_type(i1,j1)=22
     enddo

!       Consider clear condition case
!       -----------------------------
     if (ocld(1,ista)==0) then

        do ic=1,6
           if(float(abs(ocld(6+ic,ista))) < 55555) then
              write(6,*) 'cloudCover_Surface: Observed cloud above the clear level !!!'
              write(6,*) 'cloudCover_Surface: some thing is wrong in surface cloud observation !'
              write(6,*) 'cloudCover_Surface: check the station no.', ista, 'at process ', mype
              write(6,*) ic,OI(ista),OJ(ista)
              write(6,*) (ocld(k,ista),k=1,12)
              call stop2(114)
           endif
        enddo
! clean the whole column up to ceilometer height (12 kft) if ob is CLR
!    h_bk is AGL, not ASL (per Ming Hu's notes below
!
!            zlev_clr = Oelvtn(ista)+3650.
! Upcoming mods commented out below for this commit - 4/3/2010
! PH: added in column cleaning up to ceilometer height if ob is CLR
        zlev_clr = 3650.

        do k=1,nztn_p
           if (h_bk(i1,j1,k) < zlev_clr) then
              cld_cover_3d(i1,j1,k)=0.0_r_kind
              pcp_type_3d(i1,j1,k)=0
           endif
        end do

        wthr_type(i1,j1)=0

! -- Now consider non-clear obs
!    --------------------------
     else
           
!      increase zthick by 1.5x factor for ceiling < 900 m (~3000 ft - MVFR)
        cloud_dz = cloud_zthick_p
        cl_base_broken_k = -9
! ????? check with Stan    O(h_p)         if (Oelvtn(ista).lt.900.) cloud_dz = cloud_zthick_p * 2

        do ic = 1,6
           if (ocld(ic,ista)>0 .and. ocld(ic,ista)<50) then
!             if  ( csza(i,j)>=0.10 .and. sat_ctp(i1,j1)>1010.0 &
!                 .and. sat_ctp(i1,j1)<1050.)  go to 1850
!
! New tweak - 11/07/2009
!   If there was cloud in background over station but if there
!     was partial cloudiness within volume and this is one of the
!     clear columns within the polygonal area for this METAR,
!     then leave it that way and skip.
!            if (watericemax(iob,job).gt.0. .and.
!     1          kwatericemax(iob,job).gt.0 .and.
!     1          kwatericemax(iob,job).le.12) then
!              npts_cld_match = npts_cld_match + 1
!              dzbase = cl_base_ista - g3(iob,job,kwatericemax(iob,job),h_p)
!              sum_dzbase = sum_dzbase + dzbase
!              sum_dzbase_abs = sum_dzbase_abs + abs(dzbase)
!            end if

! mhu, Aug. 28, 2013: comment out patial cloudiness. It causes the degradation
! in 3000' ceiling 1-h forecast.
!              if(watericemax(ista) > 0._r_single .and. kwatericemax(i1,j1)==-1) then
!                 !PH 2/28/2013: ensure cloud building at 4 neighboring
!                 !gridpoints (Odist < 1), regardless of background 
!                 if(Odist(ista) >= 1.0_r_kind) then
!                    npts_near_clr = npts_near_clr + 1
!                    cycle   ! skip cloud build at point (i,j) because
!                            ! background is clear
!                 endif
!              endif

              if(ocld(ic,ista) == 4) then
                 if(wthr_type(i1,j1) > 10 .and. wthr_type(i1,j1) < 20) cloud_dz = 1000._r_kind  
                                         ! precipitation + highest level
                 if(wthr_type(i1,j1) == 1) cloud_dz = 10000._r_kind  ! thunderstorm
              endif

! --- calculate cloud ceiling level, not exactly, FEW SCT are also considered now
!                   iob = int(oi(ista)-idw+0.5)
!                   job = int(oj(ista)-ids+0.5)
!                   cl_base_ista = (float(ocld(6+ic,ista))+zh(iob,job))
!                   cl_base_ista = (float(ocld(6+ic,ista))+Oelvtn(ista))
! the h_bk is AGL. So observation cloud base should be AGL too, delete Oelvtn(ista)
! cover cloud base observation from AGL to ASL
              cl_base_ista = float(ocld(6+ic,ista)) + Oelvtn(ista) - zh(i1,j1)
              if(zh(i1,j1) < 1.0_r_kind .and. Oelvtn(ista) > 20.0_r_kind &
                 .and. float(ocld(6+ic,ista)) < 250.0_r_kind) then
                 cycle   ! limit the use of METAR station over oceas for low cloud base
              endif
              
              firstcloud = 0
              underlim = 10._r_kind   !

              do k=1,nztn_p
                 zdiff = cl_base_ista - h_bk(i1,j1,k)
!     Must be within cloud_dz meters (300 or 1000 currently)
!    -------------------------------------------------------------------
!  -- Bring in the clouds if model level is within 10m under cloud level.
                 if(k==1)  underlim=(h_bk(i1,j1,k+1)-h_bk(i1,j1,k))*0.5_r_kind
                 if(k==2)  underlim=10.0_r_kind    ! 100 feet
                 if(k==3)  underlim=20.0_r_kind    ! 300 feet
                 if(k==4)  underlim=15.0_r_kind    ! 500 feet
                 if(k==5)  underlim=33.0_r_kind    ! 1000 feet
                 if (k>=6 .and. k <= 7) underlim = (h_bk(i1,j1,k+1)-h_bk(i1,j1,k))*0.6_r_kind
                 if(k==8)  underlim=95.0_r_kind    ! 3000 feet
                 if(k>=9 .and. k<nztn_p-1) underlim=(h_bk(i1,j1,k+1)-h_bk(i1,j1,k))*0.8_r_kind
                 if (zdiff<underlim) then
                    if((cl_base_ista >= 1.0 .and. (firstcloud==0 .or. abs(zdiff)<cloud_dz)) .or. &
                       (cl_base_ista < 1.0 .and. (abs(zdiff)<cloud_dz)) ) then
                       if(ocld(ic,ista) == 1 ) then
                          cld_cover_3d(i1,j1,k)=max(cld_cover_3d(i1,j1,k),0.1_r_single)
                          pcp_type_3d(i1,j1,k)=0
                       elseif (ocld(ic,ista) == 2 ) then
                          cld_cover_3d(i1,j1,k)=max(cld_cover_3d(i1,j1,k),0.3_r_single)
                       elseif (ocld(ic,ista) == 3 ) then
                          cld_cover_3d(i1,j1,k)=max(cld_cover_3d(i1,j1,k),0.7_r_single)
                          if(cl_base_broken_k < 0 ) cl_base_broken_k=k
                       elseif (ocld(ic,ista) == 4 ) then
                          cld_cover_3d(i1,j1,k)=max(cld_cover_3d(i1,j1,k),1.01_r_single)
                          if(cl_base_broken_k < 0 ) cl_base_broken_k=k
                          if(wthr_type(i1,j1) == 1) then
!                            cld_type_3d(i1,j1,k)=10
                             pcp_type_3d(i1,j1,k)=1
                          endif
                          if(wthr_type(i1,j1) > 10 .and. wthr_type(i1,j1) < 20)  then
!                             cld_type_3d(i1,j1,k)=5 
                             pcp_type_3d(i1,j1,k)=1
                           endif
                       else
                           write(6,*) 'cloudCover_Surface: wrong cloud coverage observation!'
                           call stop2(114)
                       endif
                       firstcloud = firstcloud + 1
                    end if  ! zdiff < cloud_dz
                 else
!  ---- Clear up to cloud base of first cloud level
                    if (ic==1) cld_cover_3d(i1,j1,k)=0
                    if (ocld(ic,ista) == 1) pcp_type_3d(i1,j1,k)=0
                    if (ocld(ic,ista) == 3 .or. ocld(ic,ista) == 4) then
                       if( (wthr_type(i1,j1) > 10 .and. wthr_type(i1,j1) < 20)  &
                                                   .or. wthr_type(i1,j1) == 1 )  then 
                          pcp_type_3d(i1,j1,k)=1
                       endif
                    endif
                 end if  ! underlim
              end do  ! end K loop
!  ----clean cloud above stratusphere
              do k=1,nztn_p
                 if( h_bk(i1,j1,k) > 18000 ) cld_cover_3d(i1,j1,k)=0
              enddo
!
           end if     ! end if ocld > 0
        end do      ! end IC loop
!
!  clean up to broken (3) or if cloud cover less than 2, clean to cloud top
!
            if(cl_base_broken_k > 0 .and. cl_base_broken_k < nztn_p) then
                do k=1, cl_base_broken_k
                  if( cld_cover_3d(i1,j1,k) < -0.001_r_kind )   cld_cover_3d(i1,j1,k)=0
                enddo
            else
                if(ocld(1,ista) == 1 .or. ocld(1,ista) == 2 ) then
                   do k=1, nztn_p
                     if( cld_cover_3d(i1,j1,k) < -0.001_r_kind )  cld_cover_3d(i1,j1,k)=0
                   enddo
                endif
            endif

     end if      ! end if cloudy ob   ocld(1,ista) > 0

! -- Use visibility for low-level cloud whether
     if (wthr_type(i1,j1) < 30 .and. wthr_type(i1,j1) > 20 .and. &
         ocld(13,ista)  < 5000 .and. ocld(13,ista) > 1) then
           cld_type_3d(i1,j1,1) = 2
           cld_type_3d(i1,j1,2) = 2
           betav = 3.912_r_kind / (float(ocld(13,ista)) / 1000._r_kind)
           vis2qc(i1,j1) = ( (betav/144.7_r_kind) ** 1.14_r_kind) / 1000._r_kind
     endif  ! cloud or clear

   ENDDO   ! ista


!   Determine if the layer is dry or it has inversion.
!  (in either case, the cloud will be cleared out)
!
   IF(.false.) THEN     ! Set inversion strength flag
     call BckgrndCC(nlon,nlat,nsig,    &
                 t_bk,p_bk,q,h_bk,zh,  &
                 cv_bk,tbk_k,z_lcl)    ! out

     DO j = 2,nlat-1
       DO i = 2,nlon-1

         if_cloud_exist=.false.
         do k=nsig-1,2,-1
           if(cld_cover_3d(i,j,k) > 0.01_r_kind) then
             cf_model_base = cv_bk(i,j,k)
             t_model_base = tbk_k(i,j,k)
             ht_base=h_bk(i,j,k)
             if_cloud_exist=.true.
           endif
         enddo
!
! note, do we need to consider cloud base from background
         if(if_cloud_exist) then
           do k=2, nsig-1
             if(cld_cover_3d(i,j,k) > 0.01_r_kind) then
               l_cf=.false.
               l_inversion=.false.
               t_dry_adiabat = tbk_k(i,j,2) -.0098_r_kind * (h_bk(i,j,k) - h_bk(i,j,2))
               t_inversion_strength = tbk_k(i,j,k) - t_dry_adiabat
 
               IF( (tbk_k(i,j,k) > t_model_base)  .and.  &
                   (tbk_k(i,j,k) > 283.15_r_kind) .and.  &   !   temp check
                   (t_inversion_strength > 4._r_kind) ) then ! delta theta chk
                       l_inversion = .true.           ! Inversion exists
               endif
               IF( (cv_bk(i,j,k) < cf_model_base - 0.3_r_kind) .and. &
                   (h_bk(i,j,k) - ht_base >= 500._r_kind) ) THEN
                     l_cf = .true.           ! Dry layer exists
               ENDIF 
               if(l_inversion) then
                  cld_cover_3d(i,j,k) =0.0_r_kind
               endif
             endif ! in cloud
           enddo  ! k
         endif !   if_cloud_exist = true

       ENDDO   ! i
     ENDDO   ! j

   END IF     ! .true. for dry-inversion check.

END SUBROUTINE cloudCover_Surface

