!
!##################################################################
!##################################################################
!######                                                      ######
!######             SUBROUTINE INSERT_SAO1                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE insert_sao1(nx,ny,nz,dx,dy,xs,ys,zs,topo,t_sfc_k,            & 1,24
           cldcv,wtcldcv,t_modelfg,cf_modelfg,                          &
           nobsng,stn,obstype,xsng,ysng,ista_snd,                       &
           obstime,latsta,lonsta,elevsta,                               &
           kcloud,store_amt,store_hgt,                                  &
           l_stn_clouds,n_cld_snd,cld_snd,wt_snd,                       &
           i_snd,j_snd,                                                 &
           istatus)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Ingest the SAO cloud coverage observations.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  02/96.  Based on the LAPS cloud analysis code by Steve Albers,
!          1995.
!
!  MODIFICATION HISTORY:
!
!  02/06/96  J. Zhang
!            Modified for ADAS grid. Added documents.
!  03/14/97  J. Zhang
!            Clean up the code and implemented for the official
!            arps4.2.4 version
!  04/15/97  J. Zhang
!            Added error message output for the case when the
!            actual # of cloud soundings (N_CLD_SND) exceeds the
!            maximum allowed # (MAX_CLD_SND) defined in
!            adascld24.inc.
!  08/06/97  J. Zhang
!            Change adascld24.inc to adascld25.inc
!  09/09/97  J. Zhang
!            Assign a weight to a SAO cloud sounding base on
!            the cloud coverage amount.
!  09/11/97  J. Zhang
!            Change adascld25.inc to adascld26.inc
!  02/17/98  Keith Brewster
!            Changed logic to skip cloud layers with below zero
!            height -- missing, in other words.
!  05/05/98  J. Zhang
!            Abandoned the cloud grid, using the ARPS grid instead.
!  05/01/98  Keith Brewster
!            Removed abort for missing data, instead data are
!            checked for missing at every step.
!
!-----------------------------------------------------------------------
!
!  INCLUDE: (from adas.inc)
!
!  mx_sng            ! max. possible # of surface obs.
!  max_cld_snd       ! max. possible # of stations
                     ! with cloud coverage reports.
!
!  INPUT:
!
!  nx,ny,nz          ! Number of ADAS grid pts.in 3 directions.
!
!  l_stn_clouds       ! Using SAO stations' cloud obs?
!
!c ARPS grid variables.
!
!   xs     (nx)          ! The x-coord. of the physical and
                         ! computational grid. Defined at p-point.
!   ys     (ny)          ! The y-coord. of the physical and
                         ! computational grid. Defined at p-point.
!   zs    (nx,ny,nz)     ! The physical height coordinate defined at
                         ! p-point of the staggered grid.
!   topo (nx,ny)         ! The height of the terrain (equivalent
!
!   First guess fields
!
!   cf_modelfg (nx,ny,nz)  ! first guess cloud cover field
!   t_modelfg (nx,ny,nz)   ! first guess temperature field
!   t_sfc_k (nx,ny)            ! surface temperature field
!
!  cldcv (nx,ny,nz)     3D gridded fractional cloud cover analysis.
!                           (when input, it's the first guess field)
!  wtcldcv (nx,ny,nz)   weights assigned to gridded fractional
!                           cloud cover analysis. (when input, it's
!                           the weights for first guess field)
!
!   Single-level (e.g., surface) station variables
!
!   stn (mx_sng)         ! station name of single-lvel data
!   obstype (mx_sng)     ! names of sfc sources
!   obstime (mx_sng)     ! time for the observation.
!   latsta (mx_sng)      ! latitude of single-level data
!   lonsta (mx_sng)      ! longitude of single-level data
!   elevsta (mx_sng)     ! height of single-level data
!   xsng (mx_sng)        ! x location of single-level data
!   ysng (mx_sng)        ! y location of single-level data
!
!   kcloud (mx_sng)      ! number of obs. cloud layers.
!   store_amt (mx_sng,5) ! cloud coverage (ea. layer,
                         ! ea. station).
!   store_hgt (mx_sng,5) ! height of obs. cloud layers.
!
!  OUTPUT:
!
!  istatus            The flag indicating process status.
!
!  n_cld_snd         number of cloud soundings created
!  cld_snd (max_cld_snd,nz)   Cld snding obtained from SAO data.
!  wt_snd (max_cld_snd,nz)    weights assigned to cloud sounding
!                                obtained from SAO data.
!  i_snd (max_cld_snd)       i-location of each cloud sounding
!                                  station in ADAS grid
!  j_snd (max_cld_snd)       j-location of each cloud sounding
!                                  station in ADAS grid
!
!  LOCAL:
!
!  name_array (max_cld_snd)  station names (first letter) for
!                                  each cloud sounding
!  ista_snd (max_cld_snd)    index of cloud sounding stations
!  cvr_snd (max_cld_snd)     column integrated cloud cover
!  cloud_top(nx,ny)    Analyzed cloud top heights (m ASL).
!  cloud_base(nx,ny)   Analyzed cloud base heights (m ASL).
!  cloud_ceiling(nx,ny)   Analyzed cloud ceiling heights (m ASL).
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'adas.inc'
!
!-----------------------------------------------------------------------
!
!  INCLUDE: (from adas.inc)
!
!  integer mx_sng            ! max. possible # of surface obs.
!  integer max_cld_snd       ! max. possible # of stations
                             ! with cloud coverage reports.
!
!  INPUT:

  INTEGER :: nx,ny,nz

  REAL :: dx,dy              ! ADAS grid spacing
  LOGICAL :: l_stn_clouds    ! Using SAO stns' cloud obs? (or bkgrnd)

  REAL :: xs    (nx)         ! The x-coord. of the physical and
                             ! computational grid. Defined at p-point.
  REAL :: ys    (ny)         ! The y-coord. of the physical and
                             ! computational grid. Defined at p-point.
  REAL :: zs    (nx,ny,nz)   ! The physical height coordinate defined
                             ! at p-point of the staggered grid.
  REAL :: topo (nx,ny)       ! The height of the terrain (equivalent
!
!  First guess fields
!
  REAL :: cf_modelfg (nx,ny,nz)  ! first guess cloud cover field
  REAL :: t_modelfg (nx,ny,nz)   ! first guess temperature field
  REAL :: t_sfc_k (nx,ny)            ! surface temperature field
!
  REAL :: cldcv (nx,ny,nz)      ! 3D gridded frac cld cv analysis.
                                  ! (it's also bkgrnd field when input)
  REAL :: wtcldcv (nx,ny,nz)    ! wts assgned to cld cvr analysis
                                  ! (it's also bkgrnd field when input)
!
!  Surface cloud coverage reports
!
  INTEGER :: nobsng
  CHARACTER (LEN=5) :: stn (mx_sng)    ! station name of single-lvel data
  CHARACTER (LEN=8) :: obstype (mx_sng)! names of sfc sources
  INTEGER :: obstime (mx_sng)    ! time for the observation.
  REAL :: xsng (mx_sng)          ! x location of single-level data
  REAL :: ysng (mx_sng)          ! y location of single-level data
  REAL :: latsta (mx_sng)        ! latitude of single-level data
  REAL :: lonsta (mx_sng)        ! longitude of single-level data
  REAL :: elevsta (mx_sng)       ! height of single-level data
!
  INTEGER :: kcloud (mx_sng)     ! number of cloud layers.
  CHARACTER (LEN=4) :: store_amt (mx_sng,5) ! cld coverage (ea.lyr, ea. stn).
  REAL :: store_hgt (mx_sng,5)        ! height of cloud layers.
!
!  OUTPUT:
!
  INTEGER :: istatus                ! Flag for process status
!
  INTEGER :: n_cld_snd              ! # of cloud snds created
  REAL :: cld_snd (max_cld_snd,nz)   ! Cld snd obtained from SAO data.
  REAL :: wt_snd (max_cld_snd,nz)    ! wgt for SAO cld snd
  INTEGER :: i_snd (max_cld_snd)    ! i-lctn of cld snd stn in ADAS grid
  INTEGER :: j_snd (max_cld_snd)    ! j-lctn of cld snd stn in ADAS grid
!
!  LOCAL:
!
  CHARACTER (LEN=1) :: name_array (nx,ny) ! Cld snd stn name (1st letter)
  INTEGER :: ista_snd (max_cld_snd) ! index of cld snd stns
  REAL :: cvr_snd (max_cld_snd)     ! column integrated cloud cover
!
!-----------------------------------------------------------------------
!
!  Empirical factors (parameters) for cloud cover analysis
!
!-----------------------------------------------------------------------
!
!c Using SAO stations slightly outside the ARPS domain.

  INTEGER :: ix_low,ix_high,iy_low,iy_high
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  REAL :: ht_defined           ! the reachable height of the SAOs
  REAL :: ht_base,ht_top       ! base and top heights of a cloud layer
  INTEGER :: n_analyzed        ! # of stations analyzed in this routine
  INTEGER :: i,k,l
  REAL :: ri(mx_sng),rj(mx_sng)
                                ! i-,j-lctns of each stn in ADAS grid
  INTEGER :: iloc,jloc         ! i- and j-index of each sta.
  INTEGER :: igrd,jgrd       ! i- and j-index of each sta. in grid
  REAL :: cover                ! cloud fractional cover values
  REAL :: cld_thk

  INTEGER :: k_ceil

  LOGICAL :: l_out_of_bounds,l_dry,cldprt
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  istatus = 0
  WRITE(6,*)
  IF (nobsng < 1) THEN
    PRINT*,'## No SAO data available. Returning from INSERT_SAO...'
    istatus = 1
    RETURN
  END IF
  WRITE(6,*) 'Inserting SAO data from ',nobsng,' stations.'
  cldprt=(clddiag == 1)
!
  ix_low  = 1    - i_perimeter
  ix_high = nx + i_perimeter
  iy_low  = 1    - i_perimeter
  iy_high = ny + i_perimeter
!
!-----------------------------------------------------------------------
!
!  Initialize the cloud sounding array.
!
!-----------------------------------------------------------------------
!
  DO i=1,max_cld_snd
    DO k=1,nz
      cld_snd(i,k) = 0.0
      wt_snd(i,k) = 0.0
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Initially assign a unreal number for the sensor's reachable
!  height.
!
!-----------------------------------------------------------------------
!
  ht_defined = 99999.
  n_analyzed = 0
!
!-----------------------------------------------------------------------
!
!  Loop through the stations
!
!-----------------------------------------------------------------------
!
  DO i=1,nobsng

    IF( obstype(i)(1:4) == 'meso') GO TO 125
    IF( obstype(i)(1:4) == 'armmn') GO TO 125
    IF( obstype(i)(1:4) == 'isws') GO TO 125
    IF( obstype(i)(1:4) == 'coagmet') GO TO 125
    IF( obstype(i)(1:4) == 'hplains') GO TO 125
    IF( obstype(i)(1:4) == 'wpdn') GO TO 125
!
!-----------------------------------------------------------------------
!
!    Place station at proper ADAS grid point.
!
!-----------------------------------------------------------------------
!
    ri(i) = 1. + (xsng(i) - xs(1))/dx
    rj(i) = 1. + (ysng(i) - ys(1))/dy
    iloc = nint(ri(i))
    jloc = nint(rj(i))

    IF( iloc < ix_low .OR. iloc > ix_high                               &
          .OR. jloc < iy_low .OR. jloc > iy_high) THEN
!      write(6,*) 'Station is outside domain ',stn(i)
      GO TO 125
    END IF

    IF( kcloud(i) == 0) THEN        ! Kick out AMOS stations not
                                    ! reporting clouds this time
!      write(6,*)' No cloud layers - probably AMOS -',
!    :    ' goto end of loop ',stn(i),obstype(i)
      GO TO 125
    END IF

    n_analyzed = n_analyzed + 1
    n_cld_snd = n_cld_snd + 1
    IF(n_cld_snd > max_cld_snd) THEN
      PRINT*,'##ERROR: actual # of cldsnd=',n_cld_snd,                  &
             ' exceeds max. # allowed =',max_cld_snd
      PRINT*,'Please increase MAX_CLD_SND in the .inc file.'
      PRINT*,'ABORTING from INSERT_SAO......'
      STOP
    END IF

    ista_snd(n_cld_snd) = i

    IF(obstype(i)(5:5) /= ' ' .AND. obstype(i)(4:7) /= 'AMOS') THEN
!
!-----------------------------------------------------------------------
!
!      Automated Station (12000' limit)
!
!-----------------------------------------------------------------------
!
      ht_defined = elevsta(i) + 12000./3.281
    ELSE
      ht_defined = 99999.
    END IF

    IF(cldprt) THEN
      WRITE(6,1,ERR=110) stn(i),latsta(i),lonsta(i)                     &
                       ,iloc,jloc,kcloud(i)                             &
                       ,obstype(i),ht_defined
      1         FORMAT(1X,a5,' lat=',f7.2,' lon=',f7.2,' i=',i3,' j=',i3 &
               ,' kld=',i1,' obsty=',a8,' h_def=',f8.0)
      110       WRITE(6,2,ERR=3)                                        &
               (store_amt(i,l),store_hgt(i,l),l=1,kcloud(i))
      2         FORMAT(1X,5(a4,f8.0))
      3         CONTINUE
    END IF

    IF( iloc < 1 .OR. iloc > (nx-1) .OR. jloc < 1 .OR. jloc > (ny-1) ) THEN
      l_out_of_bounds = .true.
      igrd=MIN(MAX(iloc,1),(nx-1))
      jgrd=MIN(MAX(jloc,1),(ny-1))
    ELSE
      l_out_of_bounds = .false.
      igrd=iloc
      jgrd=jloc
      name_array(iloc,jloc            )=stn(i)(1:1)
      name_array(iloc,MIN(jloc+1,ny))  =stn(i)(1:1)
      name_array(iloc,MAX(jloc-1,1   ))=stn(i)(1:1)
    END IF

    cvr_snd (n_cld_snd) = 0.          ! column intg. cloud cover.
!
!-----------------------------------------------------------------------
!
!  For each station, loop through each cloud layer observed
!  by the SAO.
!
!-----------------------------------------------------------------------
!
    DO l = 1,kcloud(i)

      cover = 0.

      IF (store_hgt(i,l) < 0.0 .AND. store_amt(i,l) /= ' CLR' ) THEN
!        ht_base = cld_base(l)
!        print*,' stn',i,' level',l,' cld.ht.=',store_hgt(i,l)
!    :            ,'  adj. cld. ht.=',ht_base
!        ht_base = elevsta(i) + ht_base
        CYCLE
      ELSE
        ht_base = elevsta(i) + store_hgt(i,l)
      END IF

      IF( ht_base > ht_defined+1.) THEN

        IF( store_amt(i,l) /= ' CLR') THEN  ! Clouds

          WRITE(6,*)' Error, inconsistent SAO data, cloud base is'      &
              //' reported to be too high for this sensor type'
          WRITE(6,*) ht_base,ht_defined,obstype(i)

          WRITE(6,*)' Please check cloud layer heights in the LSO'      &
              //' file to see that they are compatable with the'
          WRITE(6,*)' types of sensors used.'
          istatus = 0
          RETURN

        ELSE                             ! CLR
          WRITE(6,*)' Note, CLR sky cloud base does not'                &
              //' reflect this sensor type'
          WRITE(6,*) ht_base,ht_defined,obstype(i)

        END IF                            ! Clouds

      END IF          !"ht_base > ht_defined+1"
!
!-----------------------------------------------------------------------
!
!  CLOUDS ARE NOW IN MSL
!  Fill in clear for entire column for SAO or up to ht_base for AWOS
!
!-----------------------------------------------------------------------
!
      IF( store_amt(i,l) == ' CLR') THEN
        cover=.01
        DO k=1,nz

          IF( zs(igrd,jgrd,k) <= ht_base                                &
                .AND. zs(igrd,jgrd,k) <= ht_defined ) THEN

            CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd           &
                         ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
          END IF
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Assign cloud cover values for entire column for SAO
!  or up to ht_base for AWOS. First check if there is dry
!  layer or inversion.
!
!-----------------------------------------------------------------------
!
      IF( store_amt(i,l) == '-SCT' ) THEN
        cover=.15 ! .2
        ht_top=ht_base+cld_thk(ht_base-elevsta(i),cover)

        DO k=2,nz-2

          IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
!            Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,0,l_dry)

            IF(.NOT. l_dry) THEN
              CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd         &
                    ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
            END IF

          ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
!            Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,1,l_dry)

          END IF  ! ht_base < zs < ht_top
        END DO    ! k = 1, nz
      END IF

      IF( store_amt(i,l) == ' SCT') THEN
        cover=.25 ! .3
        ht_top=ht_base+cld_thk(ht_base-elevsta(i),cover)

        DO k=2,nz-2

          IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
!            Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,0,l_dry)

            IF(.NOT. l_dry) THEN
              CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd         &
                    ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
            END IF

          ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
!            Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,1,l_dry)

          END IF
        END DO
      END IF

      IF( store_amt(i,l) == '-BKN' ) THEN
        cover=.4 ! .5
        ht_top=ht_base+cld_thk(ht_base-elevsta(i),cover)
        DO k=2,nz-2

          IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
!            Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,0,l_dry)

            IF(.NOT. l_dry)THEN
              CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd         &
                    ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
            END IF

          ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
!            Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,1,l_dry)

          END IF
        END DO
      END IF

      IF( store_amt(i,l) == ' BKN' ) THEN
        cover=.7
        ht_top = ht_base + cld_thk(ht_base-elevsta(i),cover) ! 1500.

        DO k=2,nz-2

          IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
!            Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,0,l_dry)

            IF(.NOT. l_dry) THEN
              CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd         &
                    ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
            END IF

          ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
!            Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,1,l_dry)

          END IF
        END DO
      END IF

      IF( store_amt(i,l) == '-OVC') THEN
        cover=.6 ! .9
        ht_top=ht_base+cld_thk(ht_base-elevsta(i),cover)
        DO k=1,nz

          IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
!            Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,0,l_dry)

            IF(.NOT. l_dry) THEN
              CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd         &
                    ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
            END IF

          ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
!            Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,1,l_dry)

          END IF
        END DO
      END IF

      IF( store_amt(i,l) == ' OVC') THEN
        cover=1.01
        ht_top = ht_base + cld_thk(ht_base-elevsta(i),cover) ! 1500.

        DO k=2,nz-2

          IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
!            Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,0,l_dry)

            IF(.NOT. l_dry)THEN
              CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd         &
                    ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
            END IF

          ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
!            Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,1,l_dry)

          END IF
        END DO
      END IF

      IF( store_amt(i,l) == ' X  ') THEN
        cover=1.01
        ht_top = ht_base + cld_thk(ht_base-elevsta(i),cover) ! 1500.

        DO k=2,nz-2

          IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
!            Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,0,l_dry)

            IF(.NOT. l_dry) THEN
              CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd         &
                    ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
            END IF
          ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
!            Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
            CALL modify_sounding(cf_modelfg,t_modelfg,topo              &
                  ,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs                      &
                  ,ht_base,ht_top,1,l_dry)

          END IF
        END DO
      END IF

      cvr_snd(n_cld_snd) = 1. - ((1. - cvr_snd(n_cld_snd)) * cover)

    END DO
!
!-----------------------------------------------------------------------
!
!    Locate the highest ceiling
!
!-----------------------------------------------------------------------

    k_ceil = nz-2
    IF( l_stn_clouds) THEN
      DO k=nz-2,2,-1
        IF( wt_snd(n_cld_snd,k) == 1.00 .AND. cld_snd(n_cld_snd,k) > 0.5) THEN
          k_ceil = k
          GO TO 1001
        END IF
      END DO
    ELSE
      DO k=nz-2,2,-1
        IF( wtcldcv(igrd,jgrd,k) == 1.00 .AND. cldcv(igrd,jgrd,k) > 0.5 ) THEN
          k_ceil = k
          GO TO 1001
        END IF
      END DO
    END IF
!
!-----------------------------------------------------------------------
!
!  Fill in other clear layers outside of clouds, below the ceiling,
!  and within defined height range of sensor.
!
!-----------------------------------------------------------------------
!
    1001    cover = .01

    IF( l_stn_clouds) THEN
      DO k=2,k_ceil
        IF( wt_snd(n_cld_snd,k) /= 1.00                                 &
              .AND. zs(igrd,jgrd,k) <= ht_defined ) THEN
          CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd             &
                ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
        END IF
      END DO
    ELSE
      DO k=2,k_ceil
        IF( wtcldcv(igrd,jgrd,k) /= 1.00                                &
              .AND. zs(igrd,jgrd,k) <= ht_defined ) THEN
          CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd             &
                ,max_cld_snd,nz,iloc,jloc,k,cover,1.)
        END IF
      END DO
    END IF

    125     CONTINUE

  END DO       ! I=1,NOBSNG
!
!-----------------------------------------------------------------------
!
  WRITE(6,*)' Num stations analyzed/cloud soundings = '                 &
                            ,n_analyzed,n_cld_snd
!
!-----------------------------------------------------------------------
!
  istatus = 1
!  999   CONTINUE
  RETURN
END SUBROUTINE insert_sao1
!
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######        SUBROUTINE MODIFY_SOUNDING                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE modify_sounding (cf_modelfg,t_modelfg,topo,                  & 14
           t_sfc_k,i_in,j_in,k,ni,nj,nk,zs,                             &
           ht_base,ht_top,init,l_dry)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Check if there is dry layer or inversion (If there is one,
!  ( the observed cloud layer will be cleared out).
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  03/96   Based on the LAPS cloud analysis code, 1995
!
!  MODIFICATION HISTORY:
!
!  02/06/96  J. Zhang
!            Modified for ADAS grid. Added documents.
!
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

!  INPUT:

  INTEGER :: init        ! flag value indicating if the layer is
                       ! within the cloud.

  INTEGER :: ni,nj               ! horizontal grid domain size
  INTEGER :: nk                  ! # of cloud analysis levels

  INTEGER :: i_in,j_in           ! i- and j- location in ADAS grid
  INTEGER :: k                   ! the cloud height level

  REAL :: cf_modelfg (ni,nj,nk)  ! first guess cloud cover field
  REAL :: t_modelfg (ni,nj,nk)   ! first guess temperature field
  REAL :: t_sfc_k (ni,nj)        ! surface temperature field
  REAL :: topo(ni,nj)            ! height of the terrain
!
  REAL :: zs(ni,nj,nk)             ! cloud analysis heights
  REAL :: ht_base,ht_top
!
!-----------------------------------------------------------------------
!
!c    OUTPUT:

  LOGICAL :: l_dry                ! if there is inversion or dry layer?
!
!-----------------------------------------------------------------------
!
!c    LOCAL:

  INTEGER :: i,j
  REAL :: cf_model_base,t_model_base,t_subcloud
  REAL :: t_dry_adiabat,t_inversion_strength

  LOGICAL :: l_wait_for_base,l_cf,l_inversion

  SAVE l_wait_for_base,cf_model_base,t_model_base,                      &
       l_inversion,t_subcloud

!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  l_dry = .false.
  l_cf = .false.
!
!-----------------------------------------------------------------------
!
!    Find ARPS grid point nearest the SAO if it is out of bounds
!
!-----------------------------------------------------------------------
!
  i = MAX(MIN(i_in,ni-1),1)
  j = MAX(MIN(j_in,nj-1),1)

  IF (init == 1) THEN ! (below base )
                      ! reset to wait for the beginning of the layer
    l_wait_for_base = .true.
    l_inversion = .false.
    t_subcloud = t_modelfg(i,j,k)
    RETURN

  ELSE ! init = 0 (inside cloud layer)
       ! For k=1 and init=0, l_wait_base=.F., and there will be NO
       ! t_subcloud available.

    IF(l_wait_for_base .OR. k <= 2) THEN
!
!-----------------------------------------------------------------------
!
!      Set reference (just within cloud base)
!
!-----------------------------------------------------------------------
!
      l_wait_for_base = .false.
      cf_model_base = cf_modelfg(i,j,k)
      t_model_base = t_modelfg(i,j,k)

!          write(6,21) t_subcloud
!21        format( /' model T_subcloud = ',f7.2)
!            write(6,1)
!1           format(1x,21x,'  cf_fg','  t_fg ',' dt_inv'
!     :                 ,' lcf',' linv','  i   j   k ',' cldht')

    END IF
!
!-----------------------------------------------------------------------
!
!  Determine if the layer is dry or it has inversion.
!  (in either case, the cloud will be cleared out)
!
!-----------------------------------------------------------------------
!
    IF(.true.) THEN     ! Set inversion strength flag
!
      t_dry_adiabat = t_sfc_k(i,j)                                      &
                   -.0098 * (zs(i,j,k) - topo(i,j))

      t_inversion_strength = t_modelfg(i,j,k) - t_dry_adiabat

      IF( ( (t_modelfg(i,j,k) > t_model_base)                           &
                            .OR.                                        &
           (k >= 2 .AND. t_modelfg(i,j,k) > t_subcloud) )               &
          .AND.                                                         &
           (t_modelfg(i,j,k) > 283.15)       & ! temp check
      .AND.                                                             &
           (t_inversion_strength > 4.)       & ! delta theta chk
      ) THEN
!
      l_inversion = .true.                  ! Inversion exists

!            write(6,2) cf_modelfg(i,j,k),t_modelfg(i,j,k)
!     :                 ,t_inversion_strength,l_cf,l_inversion
!     :                 ,i,j,k,nint(zs(i,j,k))
!
!2           format(' Inversion detected = ',3f7.2,2l4,3i4,i6)

    ELSE IF (cf_modelfg(i,j,k) < cf_model_base - 0.3   & ! cf search
        .AND. zs(i,j,k) - ht_base >= 500.) THEN
!
      l_cf = .true.           ! Dry layer exists

!                write(6,3) cf_modelfg(i,j,k),t_modelfg(i,j,k)
!     :               ,t_inversion_strength,l_cf,l_inversion
!     :               ,i,j,k,nint(zs(i,j,k))
!
!3               format(' Dry layer detected = ',3f7.2,2l4,3i4,i6)

    ELSE                          ! neither dry nor inversion
!
!-----------------------------------------------------------------------
!
!  If l_inversion and l_cf are not reset to .false. here, the WHOLE
!  cloud layer above a SINGLE level will be cleared out IF that one
!  level has inversion.
!  This is effective through the "save" statement at the beginning
!  of this subroutine.
!
!-----------------------------------------------------------------------
!
!            write(6,4) cf_modelfg(i,j,k),t_modelfg(i,j,k)
!     :              ,t_inversion_strength,l_cf,l_inversion
!     :              ,i,j,k,nint(zs(i,j,k))
!
!4           format(' model RH/T in cloud =',3f7.2,2l4,3i4,i6)

    END IF   ! inversion check

    IF( l_cf.OR.l_inversion ) THEN
      l_dry = .true.
    END IF

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

  END IF       ! INIT=1
!
!-----------------------------------------------------------------------
!
  RETURN
END SUBROUTINE modify_sounding

!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION CLD_THK                        ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION cld_thk (ht_base,cover)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Obtain thickness of a cloud layer. The thickness is simply
!  set to 1 km for cloud below 7 km and 1.5km for those above
!  7 km.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  (Jian Zhang)
!  03/96
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    ht_base          ! the height of cloud base
!    cover            ! the cloud amount of the layer
!
!  OUTPUT:
!    cld_thk          ! the thickness of the cloud layer
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: ht_base,cover,cld_thk
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF(ht_base > 7000.)THEN
    cld_thk = 1500.
  ELSE IF(ht_base > 1000.) THEN
    cld_thk = 1000.0+(ht_base-1000.)/12.
  ELSE
    cld_thk = MIN(1000.0,ht_base)
  END IF
  IF(cover > 0.01 .AND. cover < 0.5)THEN
    cld_thk = MIN(1000.,cld_thk)
  END IF

  RETURN
  END FUNCTION cld_thk

!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION CLD_BASE                       ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION cld_base (k)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Obtain the base height of a cloud layer.
!  The base height is currently an artificial function of
!  the reported cloud layer index.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  (Jian Zhang)
!  04/17/96
!
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    k          ! the index of the reported cloud layer
!
!  OUTPUT:
!    cld_base   ! the base height of the cloud layer
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  REAL :: cld_base
  INTEGER :: k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF(k == 1) cld_base = 500.
  IF(k == 2) cld_base = 2000.
  IF(k == 3) cld_base = 6000.
  IF(k == 4) cld_base = 8000.
  IF(k == 5) cld_base = 10000.
  RETURN
  END FUNCTION cld_base





!##################################################################
!##################################################################
!######                                                      ######
!######             SUBROUTINE INSERT_RADAR                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE insert_radar(nx,ny,nz,clddiag,topo,zs,temp_3d,z_lcl          & 1
           ,ref_base1,ref_base2,hgt_thresh_ref                          &
           ,grid_ref,cldcv                                              &
           ,cloud_base,cloud_base_buf,l_unresolved                      &
           ,istatus)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!    Insert radar data into cloud grid
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  05/1996
!
!  MODIFICATION HISTORY:
!
!  05/08/96  J. Zhang
!            Modified for the ARPSDAS format. Added full
!            documentation.
!  07/26/96  J. Zhang
!            Added quality check to avoid out-bounds cloud.
!  03/27/97  J. Zhang
!            Updated for the official version of arps4.2.4
!  03/28/97  J. Zhang
!            Using the lifting condensation level as the cloud base
!            when there is no SAO cloud base.
!  09/01/97  J. Zhang
!            Document cleanup for ADASCLD version 25.
!            Using the lifting condensation level as the cloud base
!            when the radar echo top is below the analyzed SAO cloud
!            base
!  09/10/97  J. Zhang
!            Put lcl in the input argument list.
!  04/20/98  J. Zhang
!            Based on the arps4.3.3 version, abandon the cloud grid.
!            Using arps grid.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
  IMPLICIT NONE

!-----------------------------------------------------------------------
!

  INTEGER :: nx,ny,nz
!  INPUT/OUTPUT:
  REAL :: cldcv(nx,ny,nz)

!  INPUT:

  INTEGER :: clddiag
  REAL :: temp_3d(nx,ny,nz)
  REAL :: zs(nx,ny,nz)
  REAL :: topo(nx,ny)
  REAL :: z_lcl(nx,ny)      ! lifting condensation level

  REAL :: grid_ref(nx,ny,nz)

  REAL :: ref_base1         ! "significant" radar echo at lower levels
  REAL :: ref_base2         ! "significant" radar echo at upper levels
  REAL :: hgt_thresh_ref    ! height criteria for "significant" radar
                            ! echo thresholds
!  OUTPUT:
  INTEGER :: istatus
  REAL :: cloud_base(nx,ny)
  REAL :: cloud_base_buf(nx,ny) ! Lowest SAO/IR base within search radius
!
!  LOCAL:
  REAL :: radar_cover
  PARAMETER(radar_cover=1.0)
  REAL :: thresh_cvr   ! lower radar echo threshold for cloud filling
  PARAMETER (thresh_cvr = 0.2)


  LOGICAL :: l_below_base
  LOGICAL :: l_unresolved(nx,ny)
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  LOGICAL :: cldprt
  REAL :: unlimited_ceiling,ref_thresh
  REAL :: zs_1d(nz)

  INTEGER :: i,j,k,kp1
  INTEGER :: icount_below,isearch_base,insert_count_tot
  INTEGER :: icount_radar_lvl,insert_count_lvl
  INTEGER :: i_out_bnd,i_lowecho
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  WRITE(6,*) 'Inserting radar reflectivity '
  cldprt=(clddiag == 1)
!
!-----------------------------------------------------------------------
!
!  Calculate Cloud Bases
!
!-----------------------------------------------------------------------
!
  unlimited_ceiling = 200000.

  DO j = 1,ny-1
    DO i = 1,nx-1
      cloud_base(i,j) = unlimited_ceiling
!
      DO k = nz-1,1,-1
        IF(cldcv(i,j,k) < thresh_cvr .AND. cldcv(i,j,k+1) >= thresh_cvr) THEN
          cloud_base(i,j) = 0.5 * (zs(i,j,k) + zs(i,j,k+1))
        END IF
      END DO ! k
!
      IF (cloud_base(i,j) > z_lcl(i,j)) THEN
        cloud_base(i,j) = z_lcl(i,j)             ! using lcl
      END IF

      cloud_base_buf(i,j) = cloud_base(i,j)
      l_unresolved(i,j) = .false.

    END DO ! i
  END DO ! j

  icount_below = 0     ! tot # of rdr echo pts below cloud base
  isearch_base = 0     ! # of neighb cloud bases being succes. found
  insert_count_tot = 0 ! tot # of data pts inserted the radar_cover
!
!-----------------------------------------------------------------------
!
!  Essentially, this go downward to detect radar tops in time
!  to search for a new cloud base
!
!-----------------------------------------------------------------------
!
  i_out_bnd=0
  i_lowecho=0

  DO i = 1,nx-1
    DO j = 1,ny-1

      DO k=1,nz
        zs_1d(k) = zs(i,j,k)
      END DO

      icount_radar_lvl = 0 ! # of lvls with refl > ref_thresh
      insert_count_lvl = 0 ! # of cloud lvls inserted radar_cover

      DO k = nz-1,1,-1
        kp1 = MIN(k+1,nz)
!
!-----------------------------------------------------------------------
!
!  Define the "significant" reflectivity value based on height (to
!  eliminate ground clutters and/or other non-weather radar echoes).
!
!-----------------------------------------------------------------------
!
        IF ((zs(i,j,k)-topo(i,j)) <= hgt_thresh_ref) THEN
          ref_thresh = ref_base1
        ELSE
          ref_thresh = ref_base2
        END IF

        IF(grid_ref(i,j,k) > ref_thresh) THEN
          icount_radar_lvl = icount_radar_lvl + 1

          l_below_base = .false.
!
!-----------------------------------------------------------------------
!
!  Test if we are at echo top
!
!-----------------------------------------------------------------------
!
          IF(k == nz-1 .OR. grid_ref(i,j,kp1) < ref_thresh) THEN
!
!-----------------------------------------------------------------------
!
!  Test if we are below the cloud base.
!
!-----------------------------------------------------------------------
!
            IF(zs(i,j,k) < cloud_base_buf(i,j)) THEN
!
!-----------------------------------------------------------------------
!
!  Radar Echo Top is below the analyzed cloud base.
!  Using the lifting condensation level as the modified cloud base
!  if it is lower than the current buffer value.
!
!-----------------------------------------------------------------------
!
              cloud_base_buf(i,j)=MIN(z_lcl(i,j),cloud_base_buf(i,j))

              IF(cloud_base_buf(i,j) < zs(i,j,k)) THEN
                isearch_base = isearch_base + 1
                IF(isearch_base < 50) THEN ! limit log output
                  WRITE(6,71) i,j,k,zs(i,j,k),cloud_base(i,j)           &
                  ,cloud_base_buf(i,j)
71                 FORMAT(' Rdr Top < Bse*gp=',3I3,' zs=',f7.0          &
                   & ,' cld_b=',f7.0,'lcl=',f7.0,' Resolved')
                END IF

              ELSE ! Probably Unresolved base
                WRITE(6,72) i,j,k,zs(i,j,k),cloud_base(i,j)             &
                ,cloud_base_buf(i,j)
72               FORMAT(1X,' **** Prob Unresolved ****'/                &
                 &  1X,'Rdr Top < Bse*gp=',3I3,' zs=',f7.0              &
                 & ,' cld_b=',f7.0,' lcl=',f7.0)

                IF(cloud_base_buf(i,j) == unlimited_ceiling) THEN
                  l_unresolved(i,j) = .true.
                  PRINT*,'Error, no LCL found for grid,',i,j,           &
                         ' aborting from INSERTRAD...'
                  STOP
                END IF ! cloud_base(i,j).eq.unlimited_ceiling
                GO TO 600

              END IF ! cloud_base_buf(i,j) .lt. zs(i,j,k)

            END IF ! Blw Cld Base: zs.lt.cloud_base_buf(i,j)

          END IF ! At Echo Top: ref(k)>thr & (k.eq.nz .or. ref(kp1)<thr
!
!-----------------------------------------------------------------------
!
!  Loop through range of cloud grid levels for this model level
!
!-----------------------------------------------------------------------
!
          IF(zs(i,j,k) > cloud_base_buf(i,j)) THEN
                             ! Insert radar if we are above cloud base
            cldcv(i,j,k) = radar_cover
            insert_count_lvl = insert_count_lvl + 1
            insert_count_tot = insert_count_tot + 1
          ELSE ! Radar Echo below cloud base
            l_below_base = .true.
          END IF

          IF(l_below_base) THEN
            icount_below = icount_below + 1

            IF(icount_below <= 50) THEN
              WRITE(6,81) i,j,k,zs(i,j,k)                               &
                          ,cloud_base_buf(i,j)
              81              FORMAT(1X,'Rdr < Bse* g.p.:',3I3          &
                                     ,' zs(i,j,k)=',f7.0,' cld_base=',f7.0)
            END IF

            GO TO 600
          END IF ! l_below_base =.true.

        ELSE ! grid_ref(i,j,k) <= ref_thresh
          IF(cldcv(i,j,k) > 0.4.AND.i_lowecho <= 25) THEN
            i_lowecho=i_lowecho+1
            IF(cldprt) WRITE(6,642) grid_ref(i,j,k),cldcv(i,j,k),i,j,k
            642            FORMAT(1X,' Rdr echo', f8.1,'<< Cld cvr',f6.2,' at',3I3)
          END IF
          GO TO 600
        END IF ! grid_ref(i,j,k) > ref_thresh ?

        IF (cldprt .AND. MOD(insert_count_tot,100) == 0) THEN
          WRITE(6,591) k,icount_radar_lvl                               &
                       ,insert_count_lvl,insert_count_tot
          591         FORMAT(1X,'Inserted radar** k=',i3/               &
                              ' n_rdrl=',i5,' n_insl=',i5,' n_instot=',i5)
        END IF
        600     CONTINUE

      END DO ! k
    END DO ! j
  END DO ! i

  WRITE(6,*)' Total cloud grid points modified by radar = '             &
             ,insert_count_tot
  istatus=1

  RETURN
END SUBROUTINE insert_radar
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE INSERT_VIS                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE insert_vis (nx,ny,nz,zs,topo,clouds_3d                       & 1
           ,albedo,cld_frac_vis_a                                       &
           ,vis_rad_thcvr,vis_rad_thdbz                                 &
           ,istat_radar,radar_ref_3d,ref_base,dbz_max_2d                &
           ,r_missing,sfc_sao_buffer,istatus)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  To process satellite data for clouds and to modify the
!  3d cloud cover array.
!  Currently Band 8 (11.2 micron) brightness temps are used with a
!  dummy call for the CO2 slicing method.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  04/26/96 (Based on LAPS code of 09/95)
!
!  MODIFICATION HISTORY:
!
!  04/26/96 (J. Zhang)
!           Modified for the ARPSDAS format.
!           Added full documentation.
!  09/02/97 (J. Zhang)
!           Added if statements to avoid putting in the reflectivity
!           threshold values in the data holes.
!  04/20/98 (J. Zhang)
!           Abandoned the cloud analysis grid, using the ARPS grid
!           instead.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz
!  INPUT/OUTPUT:
  REAL :: clouds_3d(nx,ny,nz) ! cloud cover analysis
!
!  INPUT:
  REAL :: r_missing

  REAL :: zs(nx,ny,nz)
  REAL :: topo(nx,ny)
!
!  INPUT:  parameters
  REAL :: vis_rad_thcvr          ! =0.2, defined in cloud_cv.f
  REAL :: vis_rad_thdbz          ! =10dbz, defined in cloud_cv.f
  REAL :: sfc_sao_buffer         ! =800m, defined in cloud_cv.f
!
!  INPUT: satellite visible channel data
  REAL :: albedo(nx,ny)
  REAL :: cld_frac_vis_a(nx,ny)
!
!  INPUT: radar data
  INTEGER :: istat_radar
  REAL :: dbz_max_2d(nx,ny)
  REAL :: radar_ref_3d(nx,ny,nz)
  REAL :: ref_base
!
!  OUTPUT:
  INTEGER :: istatus
!
!  LOCAL:
  INTEGER :: ih_alb(-10:20)
  INTEGER :: ih_cf_sat(-10:20)           ! Histogram for VIS cf array
  INTEGER :: ih_cfin(-10:20)             ! Histogram for input cf array
  INTEGER :: ih_cfout(-10:20)            ! Histogram for output cf array
  INTEGER :: ih_cfin_out(-10:20,-10:20)  ! Histogram for the comparison
                                       ! between input cloud cover
                                       ! array and modified cf array.
  INTEGER :: ih_cfin_sat(-10:20,-10:20)  ! Histogram for the comparison
                                       ! between input cloud cover
                                       ! array and sat. VIS cf array.
  INTEGER :: ih_cmaxin_sat(-10:20,-10:20)
  INTEGER :: ih_cmaxout_sat(-10:20,-10:20)

!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  REAL :: viscorfrc
  PARAMETER (viscorfrc=0.29)

  INTEGER :: i,j,k,kl
  INTEGER :: n_missing_albedo,n_vis_mod
  INTEGER :: i_sat,i_in,i_out,i_cmaxin,i_cmaxout
  REAL :: diffin,diffout,diffin_sum,diffout_sum                         &
       ,diffin_sumsq,diffout_sumsq
  REAL :: cld_frac_vis,cld_frac_in,cld_frac_out
  REAL :: cushion
  REAL :: cmaxin,cmaxout
  INTEGER :: iblank_radar,iset_vis
  REAL :: r_present
  LOGICAL :: l_prt

!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  istatus=0
!
!-----------------------------------------------------------------------
!
!  Initialize all histogram arrays.
!
!-----------------------------------------------------------------------
!
  DO i=-10,20
    ih_alb(i) = 0
    ih_cf_sat(i) = 0
    ih_cfin(i) = 0
    ih_cfout(i) = 0
    DO j=-10,20
      ih_cfin_out(i,j) = 0
      ih_cfin_sat(i,j) = 0
      ih_cmaxin_sat(i,j) = 0
      ih_cmaxout_sat(i,j) = 0
    END DO
  END DO
!
  n_missing_albedo = 0
  n_vis_mod = 0        ! # of data points modified by VIS data

  diffin_sum  = 0.     ! sum of diff. between column max. INPUT cld
                        ! fraction and the satellite vis cld frac'n.
  diffout_sum = 0.     ! sum of diff. between column max. OUTPUT cld
                        ! fraction and the satellite vis cld frac'n.
  diffin_sumsq  = 0.
  diffout_sumsq = 0.
!
!-----------------------------------------------------------------------
!
!  Horizontal array loop
!
!-----------------------------------------------------------------------
!
  l_prt = .false.
  DO i = 1,nx-1
    DO j = 1,ny-1
!
      kl  = nint(albedo(i,j)*10.)
      kl  = MIN(MAX(kl,-10),20)
      ih_alb(kl) = ih_alb(kl) + 1

      IF(cld_frac_vis_a(i,j) /= r_missing) THEN

        cld_frac_vis = cld_frac_vis_a(i,j)

        i_sat = nint(cld_frac_vis*10.)
        i_sat = MIN(MAX(i_sat,-10),20)
        ih_cf_sat(i_sat) = ih_cf_sat(i_sat) + 1  ! # of data points
                                  ! in a VIS cf value catagory
!
!-----------------------------------------------------------------------
!
!  Make sure satellite cloud fraction is between 0 and 1
!
!-----------------------------------------------------------------------
!
        IF(cld_frac_vis < 0.0) cld_frac_vis = 0.0
        IF(cld_frac_vis > 1.0) cld_frac_vis = 1.0

        i_sat = nint(cld_frac_vis*10.)

        cmaxin = 0.
        cmaxout = 0.

        iblank_radar = 0
        iset_vis = 0

        DO k = 1,nz-1

!          l_out(1)=k.eq.19.and.i.le.17.and.i.ge.2
!     :                    .and.j.le.43.and.j.ge.36
!          l_out(3)=k.eq.26.and.i.le.10.and.i.ge.6
!     :                    .and.j.le.26.and.j.ge.30
!          l_prt=l_out(1).or.l_out(3)

          cld_frac_in = MIN(clouds_3d(i,j,k),1.0) ! save input cf value
!
!-----------------------------------------------------------------------
!
!  Modify the cloud field with the vis input - allow .3 vis err?
!  The scheme:
!  Above the sao buffer layer, the previoud cloud analysis is
!  retained if it is within cld_frac_vis + 0.3.
!  Below the layer, the cloud analysis must be smaller than
!  cld_frac_vis. If not, reset it to cld_frac_vis.
!
!-----------------------------------------------------------------------
!
          IF(zs(i,j,k) > topo(i,j)+sfc_sao_buffer) THEN
            cushion = 0.3     ! Allow 0.3 error in VIS cld cvr.
          ELSE
            cushion = 0.0     ! VIS cld cover is the max. allowed.
          END IF

          IF(cld_frac_vis < viscorfrc .AND.                             &
                ((clouds_3d(i,j,k)-cld_frac_vis) > cushion)) THEN
            cld_frac_out = cld_frac_vis     ! VIS cf value
!
            IF(l_prt) THEN
              WRITE(6,621) i,j,k,clouds_3d(i,j,k)                       &
                           ,radar_ref_3d(i,j,k),dbz_max_2d(i,j)         &
                           ,cld_frac_vis_a(i,j)                         &
                           ,cld_frac_out
              621  FORMAT(1X,3I3,' oldcld=',f5.2,' ref=',f8.1,' dbz_m=',f8.1 &
                           ,' viscld=',f5.2,' newcld=',f5.2)
            END IF

!
!-----------------------------------------------------------------------
!
!  Determine if we need to reconcile VIS with radar
!
!-----------------------------------------------------------------------
!
            IF(istat_radar == 1 .AND. dbz_max_2d(i,j) /= r_missing      &
                  .AND. dbz_max_2d(i,j) > ref_base) THEN
                                                     ! Valid radar echo

              IF(cld_frac_out < vis_rad_thcvr) THEN

                IF(dbz_max_2d(i,j) < vis_rad_thdbz) THEN
                                  ! Blank out Radar, Normal VIS Clearing
                  iblank_radar = 1
                ELSE            ! dbz_max_2d(i,j) >= vis_rad_thdbz
                  cld_frac_out = vis_rad_thcvr     ! use radar cf
                  iset_vis = 1
                END IF

              END IF
            END IF

            IF(cld_frac_in - cld_frac_out > .01) THEN
              n_vis_mod = n_vis_mod + 1
            END IF

            clouds_3d(i,j,k) = cld_frac_out   ! Modify the output
          ELSE  ! clouds_3d - cld_frac_vis .le. cushion
            cld_frac_out = cld_frac_in        ! previous cloud analysis
          END IF  ! clouds_3d - cld_frac_vis .gt. cushion
!
!-----------------------------------------------------------------------
!
!  Update Histograms
!
!-----------------------------------------------------------------------
!
          i_in = nint(cld_frac_in*10.)
          ih_cfin(i_in) = ih_cfin(i_in) + 1

          i_out = nint(cld_frac_out*10.)
          ih_cfout(i_out) = ih_cfout(i_out) + 1

          ih_cfin_sat(i_in,i_sat)                                       &
                         = ih_cfin_sat(i_in,i_sat) + 1

          ih_cfin_out(i_in,i_out)                                       &
                         = ih_cfin_out(i_in,i_out) + 1

          cmaxin  = MAX(cmaxin,cld_frac_in)   !col.max of inp cldcvr
          cmaxout = MAX(cmaxout,cld_frac_out) !col.max of outp cldcvr

        END DO ! k
!
!-----------------------------------------------------------------------
!
!  Reconcile VIS with radar
!
!-----------------------------------------------------------------------
!
        IF(iblank_radar == 1) THEN ! NO VIS / WEAK ECHO
!
!-----------------------------------------------------------------------
!
!  Blank out radar column for this grid point
!
!-----------------------------------------------------------------------
!
          IF(cmaxout <= vis_rad_thcvr) THEN
            WRITE(6,1) i,j,cmaxout,dbz_max_2d(i,j),cld_frac_vis
            1             FORMAT(' VIS_RDR - Blank out radar: cvr/dbz/vis << ' &
                                  ,2I3,f8.2,f8.1,f8.2)
          ELSE ! Some of the column remains above the VIS threshold
               ! We are "saved" by the VIS cushion
               ! May not show up in comparisons

            WRITE(6,2) i,j,cmaxout,dbz_max_2d(i,j),cld_frac_vis
            2             FORMAT(' VIS_RDR - Blank out radar: cvr/dbz/vis-s >> ' &
                                 ,2I3,f8.2,f8.1,f8.2)
          END IF

          IF (dbz_max_2d(i,j) > ref_base)            dbz_max_2d(i,j) = ref_base
          DO kl = 1,nz
            IF (radar_ref_3d(i,j,kl) > ref_base) radar_ref_3d(i,j,kl) = ref_base
          END DO ! kl

        ELSE IF (iset_vis == 1) THEN ! NO VIS / STRONG ECHO
!
!-----------------------------------------------------------------------
!
!  Cloud cvr has been reset to threshold value above VIS
!
!-----------------------------------------------------------------------
!
          IF(cmaxout <= vis_rad_thcvr) THEN
            WRITE(6,3) i,j,cmaxout,dbz_max_2d(i,j),cld_frac_vis         &
                         ,vis_rad_thcvr
            3             FORMAT(' VIS_RDR - Reset vis: cvr/dbz/vis/thr << ' &
                                    ,2I4,f5.2,f8.1,2F5.2)

          ELSE ! Some of the column remains above the VIS threshold
               ! We are "saved" by the VIS cushion
               ! May not show up in comparisons
               ! Is resetting the VIS perhaps not necessary?

            WRITE(6,4)i,j,cmaxout,dbz_max_2d(i,j),cld_frac_vis          &
                        ,vis_rad_thcvr
            4             FORMAT(' VIS_RDR - Reset vis: cvr/dbz/vis/thr-s >> ' &
                                    ,2I4,f8.2,f8.1,2F8.2)
          END IF

        END IF ! iblank_radar=1

        IF(j == ny/2 .AND. cmaxin-cmaxout > 0.1) THEN
          WRITE(6,12) i,j,cmaxin,cmaxout
          12          FORMAT(1X,'Vismod',2I3,' colcfi=',f5.2,' colcfo=',f5.2)
        END IF

        i_cmaxin  = nint(cmaxin*10.)
        i_cmaxout = nint(cmaxout*10.)
        ih_cmaxin_sat(i_cmaxin,i_sat)                                   &
                          = ih_cmaxin_sat(i_cmaxin,i_sat) + 1
        ih_cmaxout_sat(i_cmaxout,i_sat)                                 &
                          = ih_cmaxout_sat(i_cmaxout,i_sat) + 1

        diffin  = cmaxin  - cld_frac_vis
        diffout = cmaxout - cld_frac_vis
        diffin_sum  = diffin_sum  + diffin
        diffout_sum = diffout_sum + diffout
        diffin_sumsq  = diffin_sumsq  + diffin**2
        diffout_sumsq = diffout_sumsq + diffout**2

      ELSE     ! cld_frac_vis_a(i,j) = r_missing
        n_missing_albedo =  n_missing_albedo + 1

      END IF    ! cld_frac_vis_a(i,j) .ne. r_missing

    END DO ! i
  END DO ! j

  WRITE(6,*)
  WRITE(6,*)' N_MISSING_ALBEDO = ',n_missing_albedo
    WRITE(6,*)' N_VIS_MOD = ',n_vis_mod
    WRITE(6,*)

    WRITE(6,*)'           HISTOGRAMS'
    WRITE(6,*)'   I    Alb  CFsat   CFin  CFout'
    DO i = -5,15
      WRITE(6,11)i,ih_alb(i),ih_cf_sat(i),ih_cfin(i),ih_cfout(i)
      11      FORMAT(1X,i3,4I7)
    END DO

    WRITE(6,*)
    WRITE(6,*)'  Input vs. Satellite Cloud Fraction Histogram'
    WRITE(6,*)
    WRITE(6,*)' X-axis: input cld frac, Y-axis: Satellite cld frac'
    WRITE(6,20) (i,i=0,10)
    20    FORMAT(1X,3X,11I6)
    DO j = 0,10
      WRITE(6,21) j,(ih_cfin_sat(i,j),i=0,10)
      21      FORMAT(1X,i3,11I6)
    END DO

    WRITE(6,*)
    WRITE(6,*)'  Input vs. Output Cloud Fraction Histogram'
    WRITE(6,*)
    WRITE(6,*)' X-axis: input cld frac, Y-axis: output cld fraction'
    WRITE(6,20) (i,i=0,10)
    DO j = 0,10
      WRITE(6,21) j,(ih_cfin_out(i,j),i=0,10)
    END DO

    WRITE(6,*)
    WRITE(6,*)'  Column Max Input vs. Satellite Fraction Histogram'
    WRITE(6,*)
    WRITE(6,*)' X-axis: column max.input cf, Y-axis: col.max.sat.cf.'
    WRITE(6,20) (i,i=0,10)
    DO j = 0,10
      WRITE(6,21) j,(ih_cmaxin_sat(i,j),i=0,10)
    END DO

    WRITE(6,*)
    WRITE(6,*)'  Column Max Output vs. Satellite Fraction Histogram'
    WRITE(6,*)
    WRITE(6,*)' X-axis: column max.output cf, Y-axis: col.max.sat.cf.'
    WRITE(6,20) (i,i=0,10)
    DO j = 0,10
      WRITE(6,21) j,(ih_cmaxout_sat(i,j),i=0,10)
    END DO

    r_present = nx*ny - n_missing_albedo
    IF(r_present > 0.) THEN ! write stats
      WRITE(6,31)diffin_sum/r_present,SQRT(diffin_sumsq/r_present)
      31      FORMAT(' VIS STATS: Mean/RMS input residual  = ',2F8.3)
      WRITE(6,41) diffout_sum/r_present,SQRT(diffout_sumsq/r_present)
      41      FORMAT(' VIS STATS: Mean/RMS output residual = ',2F8.3)
    END IF

    WRITE(6,*)
    istatus = 1

    RETURN
END SUBROUTINE insert_vis
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE INSERT_IR                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE insert_ir (nx,ny,nz,rlat,rlon,rland_frac,cvr_snow            & 1,5
           ,topo,zs,z_lcl,temp_3d,p_3d                                  &
           ,t_sfc_k,t_gnd_k,cldcv_sao                                   &
           ,solar_alt,solar_ha,solar_dec                                &
           ,tb8_k,tb8_cold_k,cldtop_m_tb8,cldtop_m                      &
           ,grid_spacing_m,sfc_sao_buffer                               &
           ,cldcv,istatus)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  To process satellite data for clouds and to modify the
!  3d cloud cover array.
!  Currently Band 8 (11.2 micron) brightness temps are used.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  (Jian Zhang)
!  04/1996  Based on the LAPS cloud analysis code (Steve Albers,
!           09/1995)
!
!  MODIFICATION HISTORY:
!
!  04/26/96  J. Zhang
!            Modified for the ARPSDAS format.
!            Added full documentation.
!  07/19/96  J. Zhang
!            Added the quality check part to deal with the tall
!            and steep cloud towers.
!  07/24/96  J. Zhang
!            Added the quality check for cases when the cloud top
!            is higher than the model top boundary
!  11/01/96  J. Zhang
!            Added the rawinsonde data for determing the cloud
!            top height for the low clouds.
!  03/18/97  J. Zhang
!            Cleanup the code and implemented for the official
!            arps4.2.4 version
!  08/06/97  J. Zhang
!            Change adascld24.inc to adascld25.inc
!  09/09/97  J. Zhang
!            Fixed a bug when calling COMPARE_RAD. Further cleanup.
!  09/10/97  J. Zhang
!            Lifting condensation level is now an input argument.
!            Change adascld25.inc to adascld26.inc
!  04/20/98  J. Zhang
!            Abandoned the cloud analysis grid, using the ARPS grid
!            instead.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:

  INCLUDE 'phycst.inc'
  INCLUDE 'adas.inc'
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz            ! ARPS grid size
!  INCLUDE: ( from adas.inc)
!  real r_missing
!
!  INPUT/OUTPUT:
  REAL :: cldcv(nx,ny,nz)        ! 3D cldcv analx after useing sat.data
!
!  INPUT:   general
  REAL :: grid_spacing_m
!
!  INPUT:  for model grid geometry
  REAL :: topo(nx,ny)
  REAL :: rlat(nx,ny),rlon(nx,ny) ! Lat. and lon. of the grid pt.
  REAL :: zs(nx,ny,nz)            ! Phys. ht (m) at model grid point.
  REAL :: z_lcl(nx,ny)            ! Lifting condensation level (MSL)
  REAL :: temp_3d(nx,ny,nz)       ! temp. analysis on model grid
  REAL :: p_3d(nx,ny,nz)          ! pres. analysis on model grid
!
!  INPUT: for model grid sfc characteristic fields
  REAL :: rland_frac(nx,ny)  ! land coverage fraction at a g.p.
  REAL :: cvr_snow(nx,ny)         ! snow coverage at a grid point.
  REAL :: solar_alt(nx,ny)        ! solar altitude angle
  REAL :: solar_ha(nx,ny)         ! solar hour angle
  REAL :: solar_dec
!
!  INPUT: for model grid analysis fields
  REAL :: t_sfc_k(nx,ny)          ! sfc air temp.
!
!  INPUT:  for SAO
  REAL :: sfc_sao_buffer          ! No clearing cloud by sat. below
                                  ! this ht.(m AGL) (hence letting
                                  ! SAOs dominate)
  REAL :: cldcv_sao(nx,ny,nz) ! 3D Cloud cover array
!
!  OUTPUT:
  INTEGER :: istatus
  REAL :: t_gnd_k(nx,ny)          ! ground skin temp.
  REAL :: tb8_k(nx,ny)            ! Obs. sate. band8 bright. temp.
  REAL :: tb8_cold_k(nx,ny)       ! Cold filtered band8 bright. temp.
  REAL :: cldtop_m(nx,ny)         ! Adj. Sat. cloud top height (m)
  REAL :: cldtop_m_tb8(nx,ny)     ! Sat. cld top ht (m) from band 8 data
!
!  Local: factors and parameters
  REAL :: sfc_ir_buffer           !No adding cld by sat. below this ht.
                                  ! (m AGL) (hence letting SAOs dominate)
  PARAMETER (sfc_ir_buffer = 3000.)

  REAL :: thk_def                 !Def. cld thkness inserted by sat.data
  PARAMETER (thk_def = 1500.)

  REAL :: thr_sao_cvr             ! cldcvr thres. in evaluating presence
                                  ! of SAO cld layers
  PARAMETER (thr_sao_cvr = 0.1)

  REAL :: thresh2                 ! Thres.for IR cloud detection
  PARAMETER (thresh2 = 3.5)    ! (Tsfc - T_IR) Were 5K and 15K also
!
!  LOCAL: for band 8 brightness temp.
  REAL :: tb8_est                 ! estimated tb8 temp. by using
                                  ! cld cvr analysis
!  LOCAL: for interpolations
  REAL :: zs_1d(nz)
  REAL :: t_1d(nz)  ! 3D temperature analysis on model grid
  REAL :: p_1d(nz)  ! pres. analysis on model grid
!
!  LOCAL: for Controling search box for SAO analyzed data
  INTEGER :: idelt(3)
  INTEGER :: jdelt(3)
!
!  LOCAL: for satellite cloud presence and cloud top height
  LOGICAL :: l_tb8,l_cloud_present
  REAL :: cldtop_m_old
  INTEGER :: nlyr               ! for iterative adj. of cld cvr field
!
!  FUNCTIONS:
  REAL :: t_ground_k
  REAL :: temp_to_rad
  REAL :: rad_to_temp
  REAL :: band8_cover
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,i_sao,j_sao
  INTEGER :: il,ih,jl,jh,ii,jj
  INTEGER :: n_no_sao1,n_no_sao2,n_no_sao3
  INTEGER :: i_reset_base

  REAL :: htbase,htbase_init,ht_sao_base,ht_cloud
  PARAMETER (ht_cloud = 2000.0)

  INTEGER :: i_delt,idelt_max,idelt_index,jdelt_index

  REAL :: sat_cover,temp_thresh,cldcv_above                             &
       ,cover,cover_new,buffer
  REAL :: tmin
  INTEGER*4 nidelt,njdelt
  DATA nidelt/3/,njdelt/3/

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Print some band 8 brightness temp. samples.
!
!-----------------------------------------------------------------------
!
  WRITE(6,*) 'Insert satellite IR data routines called.'
  WRITE(6,*)
!
  il=MAX(1,nx/2-5)
  ih=MIN(nx,nx/2+5)
  jl=MAX(1,ny/2-5)
  jh=MIN(ny,ny/2+5)

  WRITE(6,*) ' Band 8 brightness temp values in the mid-domain:'
  WRITE(6,601) (i,i=il,ih)
  601   FORMAT(/1X,3X,11(2X,i2,2X))
  WRITE(6,603) (j,(tb8_k(i,j),i=il,ih),j=jh,jl,-1)
  603   FORMAT(1X,i3,11F6.1)
!
!-----------------------------------------------------------------------
!
!  Define a box for searching SAO cloud base
!
!-----------------------------------------------------------------------
!
  idelt_max = nint(50000. / grid_spacing_m)

  idelt(1) = -idelt_max
  idelt(2) = 0
  idelt(3) = +idelt_max

  jdelt(1) = 0
  jdelt(2) = +idelt_max
  jdelt(3) = -idelt_max
!
!-----------------------------------------------------------------------
!
!  First guess conversion from cloud height grid to ADAS grid
!  This has to err slightly on the high side
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!  QC check on band 8 (11.2 mm) brightness temps
!
!-----------------------------------------------------------------------
!
  DO j = 1,ny
    DO i = 1,nx
      IF( tb8_k(i,j) < 173. .OR. tb8_k(i,j) > 350.) tb8_k(i,j)=-999.
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate cold filtered temperatures
!
!-----------------------------------------------------------------------
!
  i_delt = MAX(1,nint(grid_spacing_m/10000.))
  DO j = 1,ny
    jl = MAX(1 ,j-i_delt)
    jh = MIN(ny,j+i_delt)

    DO i = 1,nx
      il = MAX(1 ,i-i_delt)
      ih = MIN(nx,i+i_delt)
      tb8_cold_k(i,j) = 1.0E15

      DO jj = jl,jh,i_delt
        DO ii = il,ih,i_delt
          IF (tb8_k(ii,jj) > 0.)                                        &
              tb8_cold_k(i,j) = MIN(tb8_cold_k(i,j),tb8_k(ii,jj))
        END DO ! ii
      END DO ! jj

      IF(tb8_cold_k(i,j) > 350.) tb8_cold_k(i,j) = -999.

    END DO ! i
  END DO ! j
!
!-----------------------------------------------------------------------
!
!  Calculate ground temperature (for now equate to sfc air temp)
!
!-----------------------------------------------------------------------
!
  WRITE(6,*) 'Computing ground skin temperature'
  DO j = 1,ny-1
    DO i = 1,nx-1
      IF ( rland_frac(i,j) >= 0.5) THEN
        t_gnd_k(i,j) = t_ground_k(t_sfc_k(i,j),solar_alt(i,j)           &
                                  ,solar_ha(i,j),solar_dec              &
                                  ,rlat(i,j),cvr_snow(i,j)              &
                                  ,r_missing,i,j,nx,ny)
      ELSE ! water environment
        t_gnd_k(i,j) = t_sfc_k(i,j)
      END IF
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Check number of points thrown out as a function of offset to check
!  for IR navigation errors. This is only a diagnostic process,
!  does not change anything.
!
!-----------------------------------------------------------------------
!
  CALL correlation(t_gnd_k,tb8_k,thresh2,nx,ny)
!
!-----------------------------------------------------------------------
!
!  Initializes routine
!
!-----------------------------------------------------------------------
!
  n_no_sao1 = 0 ! # of g.p with sat.cld but no SAO cld at the grid pt
  n_no_sao2 = 0 ! # of g.p with sat.cld but no SAO cld in the box
  n_no_sao3 = 0 ! # of pts with sat.cldtop lower than SAO cld base.
  i_reset_base = 1  ! # of pts that sate. ceiling being resetted
!
!-----------------------------------------------------------------------
!
  WRITE(6,*) 'Find cloud_top for each grid point.'

  DO j=1,ny-1
    DO i=1,nx-1
      IF(tb8_k(i,j) > 0.) THEN
        DO k=1,nz
          zs_1d(k) = zs(i,j,k)
          t_1d(k)  = temp_3d(i,j,k)
          p_1d(k)  = p_3d(i,j,k)
        END DO
        tmin=temp_3d(i,j,2)
        DO k=2,nz-1
          tmin=MIN(tmin,temp_3d(i,j,k))
        END DO
!
!-----------------------------------------------------------------------
!
!  Calculate cloud top height from Band 8 method
!
!-----------------------------------------------------------------------
!
        CALL cloud_top(nx,ny,nz,                                        &
                    i,j,zs_1d,t_1d,p_1d,z_ref_lcl,z_lcl(i,j),           &
                    topo(i,j),tmin,r_missing,                           &
                    tb8_k(i,j),t_gnd_k(i,j),thresh2,                    &
                    cldtop_m_tb8(i,j),cldtop_m(i,j),                    &
                    l_tb8,l_cloud_present,sat_cover)
!
!-----------------------------------------------------------------------
!
!  Modify those levels where the satellite shows warmer than the
!  calculated brightness temp from the analysis of SAO/Pireps
!
!-----------------------------------------------------------------------
!
        DO k = nz-1,1,-1

          IF (cldcv(i,j,k) > 0.04) THEN ! Efficiency test
!
!-----------------------------------------------------------------------
!
!  Estimate the brightness temperature using the cloud cover
!  field.
!
!-----------------------------------------------------------------------
!
            IF (cldcv(i,j,k) /= r_missing) THEN
              tb8_est = temp_to_rad(temp_3d(i,j,k))                     &
                  * cldcv(i,j,k) + temp_to_rad(t_gnd_k(i,j))            &
                  * (1.-cldcv(i,j,k))
              tb8_est = rad_to_temp(tb8_est)  ! estimated tb8 temp.

!                tb8_est_a = temp_3d(i,j,k) * cldcv(i,j,k)
!     :             + t_gnd_k(i,j) * (1.-cldcv(i,j,k))

            ELSE
              tb8_est = t_gnd_k(i,j)
            END IF
!
!-----------------------------------------------------------------------
!
!  Test if clouds detected by SAO/PIREP should have been
!  detected by the satellite (if satellite is warmer than analysis)
!
!-----------------------------------------------------------------------
!
            IF (tb8_k(i,j) > tb8_est) THEN !seems too much SAO cld

!         Don't touch points within buffer of surface
              IF (zs(i,j,k)-topo(i,j) > sfc_sao_buffer) THEN

!         Does satellite still imply at least some cloud?
                IF(t_gnd_k(i,j)-tb8_k(i,j) > thresh2) THEN ! Some cloud
                  IF(cldcv(i,j,k) > 0.9) THEN
                    cldcv(i,j,k)=.01           ! Lower top of solid cld
                  ELSE                        ! Cover < 0.9, correct it
                    cldcv(i,j,k) = band8_cover( tb8_k(i,j)              &
                        ,t_gnd_k(i,j),temp_3d(i,j,k))

                    cldcv(i,j,k) = MAX(0.0,MIN(1.0,cldcv(i,j,k)))
                  END IF

                ELSE !tg-tb<8 Band 8 nearly matches grnd, clear it
!
!-----------------------------------------------------------------------
!
!           Insure that "black (or grey) stratus" is not present
!
!-----------------------------------------------------------------------
!
                  temp_thresh = MIN(t_gnd_k(i,j),t_sfc_k(i,j)-10.0)
                  IF(temp_3d(i,j,k) < temp_thresh)THEN
                    cldcv(i,j,k)=.01 ! not in inversion, clear it out
                  END IF

                END IF ! tg-tb>8 IR signature present
              END IF ! cldht-topo.gt.sfc_sao_buffer is .true.
            END IF ! tb8_k(i,j).gt.tb8_est
          END IF ! Current Cloud Cover is significant (> .04)
        END DO ! k (for clearing clouds)
!
!-----------------------------------------------------------------------
!
!  Insert satellite clouds
!
!-----------------------------------------------------------------------
!
        IF(l_cloud_present) THEN ! Insert satellite clouds
!
!-----------------------------------------------------------------------
!
!  Set initial satellite cloud base.
!
!-----------------------------------------------------------------------
!
          htbase_init=cldtop_m(i,j) - thk_def
          htbase = htbase_init
!
!-----------------------------------------------------------------------
!
!  Locate lowest SAO cloud base
!
!-----------------------------------------------------------------------
!
          ht_sao_base = 1E30
          cldcv_above = cldcv_sao(i,j,nz-1)

          DO k=nz-2,1,-1

            IF (cldcv_sao(i,j,k) <= thr_sao_cvr                         &
                  .AND. cldcv_above > thr_sao_cvr) THEN
              ht_sao_base = zs(i,j,k+1)
            END IF

            cldcv_above = cldcv_sao(i,j,k)
          END DO ! k

          i_sao = i
          j_sao = j

!jz         IF (ht_sao_base.eq.1e30) THEN ! Srch for nearby SAO cld lyrs
          IF (ht_sao_base >= 1E20) THEN ! Srch for nearby SAO cld lyrs
                                        ! because the sat.says cld
                                        ! and the SAO doesn't
            n_no_sao1 = n_no_sao1 + 1 ! Sat.cld but no SAO cld at same pt.

            DO jdelt_index = 1,njdelt
              jj = j + jdelt(jdelt_index)  ! jdelt: 0,5,-5

              DO idelt_index = 1,nidelt
                ii = i + idelt(idelt_index)  ! idelt: -5,0,5

                IF(ii >= 1.AND.ii <= nx-1 .AND. jj >= 1.AND.jj <= ny-1) THEN
                  cldcv_above = cldcv_sao(ii,jj,nz-1)

                  DO k = nz-2,1,-1

                    IF (cldcv_sao(ii,jj,k) <= thr_sao_cvr               &
                          .AND. cldcv_above > thr_sao_cvr) THEN
                      ht_sao_base = zs(i,j,k+1)
                    END IF

                    cldcv_above = cldcv_sao(ii,jj,k)
                  END DO ! k

!jz               if(ht_sao_base .ne. 1e30)then
                  IF(ht_sao_base < 1E20)THEN
                    i_sao = ii
                    j_sao = jj
                    goto 201
                  END IF

                END IF  ! In bounds
              END DO ! IDELT_INDEX
            END DO ! JDELT_INDEX
          END IF  ! ht_sao_base.eq.1e30

!jz201      IF (HT_SAO_BASE .EQ. 1E30) THEN ! Sat. cld but no SAO cld

201         IF (ht_sao_base >= 1E20) THEN ! Sat. cld but no SAO cld
            n_no_sao2 = n_no_sao2 + 1     ! even in neighbor points.
            cover=sat_cover
            htbase_init = ht_sao_base

            IF(tb8_k(i,j)-t_gnd_k(i,j) < -21.) THEN ! We more likely
                                                    ! have a cloud
              buffer = 2100.
            ELSE
              buffer = sfc_ir_buffer ! Weed out IR tops < ~5000m AGL
            END IF
!
!-----------------------------------------------------------------------
!
!  Calculate new cloud top and cover
!
!-----------------------------------------------------------------------
!
            cldtop_m_old = cldtop_m(i,j)
!
            CALL cloud_top(nx,ny,nz,                                    &
                    i,j,zs_1d,t_1d,p_1d,z_ref_lcl,z_lcl(i,j),           &
                    topo(i,j),tmin,r_missing,                           &
                    tb8_cold_k(i,j),t_gnd_k(i,j),thresh2,               &
                    cldtop_m_tb8(i,j),cldtop_m(i,j),                    &
                    l_tb8,l_cloud_present,sat_cover)
!
!-----------------------------------------------------------------------
!
!  Change to cover
!
!-----------------------------------------------------------------------
!
            cover=band8_cover(tb8_k(i,j),t_gnd_k(i,j),tb8_cold_k(i,j))
            htbase = MAX(topo(i,j)+buffer, cldtop_m(i,j)-thk_def )
            i_reset_base = i_reset_base +1

            IF (htbase > cldtop_m(i,j)) THEN
              n_no_sao3 = n_no_sao3 + 1 ! Sat.cld, but no SAO cld,
                                       ! & sat.cld is too low.
            ELSE IF(i_reset_base/50*50 == i_reset_base) THEN
              WRITE(6,611) i,j
              611             FORMAT(/1X,' Correction at point(',2I2    &
                                   ,') because there is low sat.cld,'   &
                                   ,' but no SAO cld.')
              WRITE(6,612,ERR=212) tb8_k(i,j),tb8_cold_k(i,j)           &
                   ,cldtop_m_old,cldtop_m(i,j),cover
              612             FORMAT(1X,'tb=',f8.1,'  tb_cold=',f8.1,   &
                                    ' ctop=',f9.0,' ctop_cold=',f9.0,' cvr=',f5.2/)
            212           END IF

          ELSE IF (ht_sao_base > cldtop_m(i,j)) THEN ! Sat.cld, nut no
                                     ! SAO cld at same pt. Do have
                                     ! SAO cld in nearby pt. AND
                                     ! Sat.top is below ceiling
            cover=sat_cover
            htbase_init = ht_sao_base   ! using SAO cld base
            htbase = htbase_init
            cldtop_m_old = cldtop_m(i,j)
            cldtop_m(i,j) = htbase_init + thk_def ! correct sat. cldtop
!
!-----------------------------------------------------------------------
!
!  Find a thinner value for cloud cover consistent with the new
!  higher cloud top and the known brightness temperature.
!  Note that cover is not really used here as an input
!
!-----------------------------------------------------------------------
!
            CALL correct_cover(cover,cover_new                          &
                    ,cldtop_m_old,cldtop_m(i,j),i,j,nx,ny,nz            &
                    ,zs_1d,t_1d,tb8_k(i,j),t_gnd_k(i,j)                 &
                    ,istatus)
!
            IF (istatus /= 1) THEN
              WRITE(6,*)' Error in correct_cover'
              RETURN
            END IF
            cover = cover_new

          ELSE ! Normal use of satellite data
               ! There is a ht_sao_base below cldtop_m in the area
               ! near the grid pt.

            cover=sat_cover
!
!-----------------------------------------------------------------------
!
!  Locate SAO cloud base below satellite cloud top, modify
!  satellite cloud base. Highest SAO ceiling within default
!  thickness range of satellite layer is used.
!
!-----------------------------------------------------------------------
!
            DO k=nz-1,1,-1

              IF (zs(i,j,k) >= htbase_init .AND.                        &
                    zs(i,j,k) <= cldtop_m(i,j)) THEN

                IF (cldcv_sao(i_sao,j_sao,k) <= thr_sao_cvr .AND.       &
                      cldcv_sao(i_sao,j_sao,k+1) > thr_sao_cvr) THEN

!c            We have an SAO base

                  htbase = zs(i,j,k+1)

!c              If SAO (hence satellite) base is above the satellite
!c              cloud top, lower the satellite base by one grid level

                  IF (htbase > cldtop_m(i,j)) htbase = zs(i,j,k)

                  GO TO 301

                END IF ! We have an SAO base
              END IF ! in satellite layer
            END DO ! k
          END IF ! HT_SAO_BASE .EQ. 1E30
                 ! No SAO cloud base in the area near the grid pt.

          301         IF (htbase /= htbase_init  .AND.                  &
                            i_reset_base/50*50 == i_reset_base) THEN
            WRITE(6,*)' Satellite ceiling reset by SAO.'
            WRITE(6,660) i,j,htbase,htbase_init,cldtop_m(i,j)
            660           FORMAT(1X,'i/j=',2I3,' htbase_new/old=',2F7.0, &
                               ' cldtop_ir=',f7.0)
          END IF
!
!-----------------------------------------------------------------------
!
!  Add satellite cloud to array
!
!-----------------------------------------------------------------------
!
          DO k=nz,1,-1
            IF (zs(i,j,k) >= htbase .AND.                               &
                zs(i,j,k) <= cldtop_m(i,j)) & ! in satellite layer
            cldcv(i,j,k)=cover
          END DO

        END IF ! l_cloud_present (Cloudy)

      END IF ! non-missing IR data

    END DO ! i=1,nx
  END DO ! j=1,ny
!
  istatus = 1
!
!-----------------------------------------------------------------------
!
!  Write stats on Band 8 (11.2mm) methods
!
!-----------------------------------------------------------------------
!
  WRITE(6,*)' n_no_sao (1/2/3) = ',n_no_sao1,n_no_sao2,n_no_sao3

  CALL compare_rad (nx,ny,nz,clddiag,r_missing,zs,temp_3d               &
        ,cldcv,t_sfc_k,t_gnd_k,tb8_k                                    &
        ,cvr_snow,nlyr)
!
!-----------------------------------------------------------------------

!  999   RETURN
END SUBROUTINE insert_ir





!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE CORRELATION                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE correlation(t,tb8_k,thresh,ni,nj) 1
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Checking the navigation of the IR satellite data
!  It's only a diagnostic procedure.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Dan Birkenheuer?)
!  09/95.
!
!  MODIFICATION HISTORY:
!
!  04/29/96 J. Zhang
!           Modified for the ARPSDAS format.
!           Added full documentation.
!
!  05/01/98 Keith Brewster
!           Added check for missing data.
!
!-----------------------------------------------------------------------
!
!  Variable Declaration
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!  INPUT:
  INTEGER :: ni,nj      ! Model grid size
  REAL :: t(ni,nj)      ! ground skin temperature
  REAL :: tb8_k(ni,nj)  ! satellite IR band 8 brightness temperature
  REAL :: thresh        ! threshold define the bias criteria
!
!  LOCAL:
  INTEGER :: ibox
  PARAMETER (ibox = 3)
  INTEGER*4 i_corr_array(-ibox:ibox,-ibox:ibox)
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: min_count,ioff,joff,i,j,ii,jj,icount                       &
      ,ioff_min,joff_min
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  WRITE(6,*)
  WRITE(6,*)' Checking Navigation of IR satellite, thresh = '           &
                                            ,thresh

  min_count = 999999

  DO joff = -ibox,ibox
    DO ioff = -ibox,ibox
      icount = 0

      DO j = 1,nj-1
        jj = MAX(MIN(j+joff,nj-1),1)

        DO i = 1,ni-1
          ii = MAX(MIN(i+ioff,ni-1),1)

          IF(tb8_k(ii,jj) > 0. .AND. (t(i,j) - tb8_k(ii,jj)) > thresh)  &
              icount = icount+1 ! Clouds Indicated

        END DO ! i
      END DO ! j

      IF(icount < min_count)THEN
        min_count = icount
        ioff_min = ioff
        joff_min = joff
      END IF

      i_corr_array(ioff,joff) = icount

    END DO ! IOFF
  END DO ! JOFF

  WRITE(6,99)
  99    FORMAT(//1X,' Number of cloudy points indicated'                &
               ,' by offseting tb8:')
  WRITE(6,100) (i,i=-ibox,ibox)
  100     FORMAT(1X,'joff:ioff',7I5)
  DO j = ibox,-ibox,-1
    WRITE(6,101) j,(i_corr_array(i,j),i=-ibox,ibox)
    101     FORMAT(1X,1X,i2,5X,7I5)
  END DO

  WRITE(6,*)
  WRITE(6,201) min_count,ioff_min,joff_min
  201   FORMAT(' Min. of',i5,' pts flagged with an offset of',2I4/)

  RETURN
END SUBROUTINE correlation




!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE CLOUD_TOP                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cloud_top(nx,ny,nz,                                          & 2
           i,j,zs_1d,t_1d,p_1d,z_ref_lcl,z_lcl,                         &
           topo,tmin,r_missing,                                         &
           tb8_k,t_gnd_k,thresh2,                                       &
           cldtop_m_tb8,cldtop_m,                                       &
           l_tb8,l_cloud_present,sat_cover)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This routine computes the cloud top height given a band 8
!  brightness temperature and 3D fields of temp and height.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  04/96 Based on the LAPS cloud analysis code (09/95 Dan Birkenheuer)
!
!  MODIFICATION HISTORY:
!
!  04/26/96 (J. Zhang)
!  Modified for the ARPSDAS format.
!  Added full documentation.
!
!  11/01/96 (J. Zhang)
!  Added a scheme for calculation of cloud top height of low clouds
!  (e.g., persistent stratocumulus), where temperature inversion
!  may exist and the simple matching scheme may fail.
!  Reference:
!  Macpherson et al., 1996: The impact of MOPS moisture data in
!      the U.K. Meteorological Office mesoscale data assimilation
!      scheme.   MWR, vol. 124,  1746-1766
!
!  09/10/97 (J. Zhang)
!           Lifting condensation level is input through calling
!           argument.
!
!-----------------------------------------------------------------------
!
!  Variable Declarition
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  For tb8 method
!
!-----------------------------------------------------------------------
!
!  INPUT:
  INTEGER :: nx,ny,nz     ! The model grid size
!
!  INPUT: Vertical 1D arrays
  REAL :: zs_1d(nz)    ! physical height (m) at scalar points
  REAL :: t_1d(nz)  ! 3D temperature analysis on model grid
  REAL :: p_1d(nz)  ! pres. analysis on model grid
!
!  INPUT: At grid point (i,j)
  INTEGER :: i,j          ! grid point location indices
  REAL :: z_lcl           ! Lifting condensation level at i,j
  REAL :: tb8_k           ! band 8 beightness temperature(K)
  REAL :: t_gnd_k         ! gound temperature
  REAL :: topo            ! terrain height at the grid point
!
!  INPUT: Constants
  REAL :: tmin
  REAL :: r_missing       ! missing data flag value
  REAL :: thresh2         ! Input from parent routine (=8K)
  REAL :: z_ref_lcl       ! ref. level for computing LCL
!
!  OUTPUT:
  LOGICAL :: l_cloud_present ! "cloudy" indicated by sate. IR data?
  LOGICAL :: l_tb8           ! "cloudy" indicated by band 8 data?
  REAL :: cldtop_m_tb8       ! cld top height(m) obtained from tb8_k
  REAL :: cldtop_m           ! cld top height(m) obtained from IR data
  REAL :: sat_cover          ! cloud fraction obtained from IR data
!
!  LOCAL:
  REAL :: gamma_s         ! ref. moist adiabatic lapse rate
!
!  LOCAL:  For the conceptual model for the cloud top
  REAL :: t_base,p_base_pa,zbase,ztop,dzinc,zht
  PARAMETER(dzinc=100.0)
  REAL :: p,press,temp,eso,dlvdt,des,dtz,rl,es
  REAL :: tb8_loc
  INTEGER :: nlevel,nlm1
!
!  CONSTANTS:
  REAL :: gamma_d   ! dry adiabatic lapse rate (K/m)
  PARAMETER (gamma_d = 9.8/1004.0)
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: kl,j1
  REAL :: arg,frac_k,frac_z,z_ref,t_ref
!
!-----------------------------------------------------------------------
!
!  Include files.
!
  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Define constants
!
!-----------------------------------------------------------------------
!
  dlvdt=-2.3693E+3        ! J/kg/K
  eso=610.78              ! pa
!
!-----------------------------------------------------------------------
!
!  This section finds the cloud top using Band 8 data
!
!-----------------------------------------------------------------------
!
  cldtop_m_tb8 = r_missing ! or zeros
  l_tb8 = .false.
!
  IF (tb8_k - t_gnd_k < -thresh2) THEN ! probably clds
    l_tb8 = .true.
    IF (tb8_k < 253.15) GO TO 300 ! Using scan&match method
!
!-----------------------------------------------------------------------
!
!  When Tb >= -20 deg C, find the cloud top height using a
!  conceptual model presented in Macpherson et al. (1996).
!
!-----------------------------------------------------------------------
!
    zbase = z_lcl
    IF (zbase >= zs_1d(nz-1)) GO TO 300
! Using scanning & matching method
!
!-----------------------------------------------------------------------
!
!  Find the model pressure and temperature at the lifting
!  condensation level.
!
!-----------------------------------------------------------------------
!
    DO kl = 2,nz-1
      z_ref = z_ref_lcl + topo
      IF (z_ref <= zs_1d(kl)) THEN
        frac_z = (z_ref-zs_1d(kl-1))/(zs_1d(kl)-zs_1d(kl-1))
        t_ref = t_1d(kl-1) + frac_z * (t_1d(kl) - t_1d(kl-1))
        GO TO 50
      END IF
    END DO  ! kl
    PRINT*,'Error for Tbase,aborting....'
    STOP
    50      CONTINUE
    t_base = t_ref - (z_lcl -z_ref)*gamma_d
!
    DO kl = 2,nz-1
      IF (zbase <= zs_1d(kl)) THEN
        frac_z = (zbase-zs_1d(kl-1))/(zs_1d(kl)-zs_1d(kl-1))
        arg = ALOG(p_1d(kl-1))                                          &
               + frac_z * (ALOG(p_1d(kl))-ALOG(p_1d(kl-1)))
        p_base_pa = EXP(arg)
        GO TO 100
      END IF
    END DO  ! kl
    PRINT*,'Error for Zbase,aborting....'
    STOP
    100     CONTINUE
    IF (t_base <= tb8_k) GO TO 300
! Using scanning & matching method
!
!-----------------------------------------------------------------------
!
!  Find the cloud top height, which is the level where the temp.
!  reaches the satellite brightness temperature throught the moist
!  adiabatic cooling process.
!
!-----------------------------------------------------------------------
!
    temp   = t_base    !  temp   = t_base_k
    press  = p_base_pa
    gamma_s = 0.004        ! K/m
    nlevel = (temp -tb8_k)/gamma_s/dzinc  ! a guess for cloud top
    IF(nlevel < 0) THEN
      PRINT*,'Error in zlevel. aborting....'
      PRINT*,' nlevel=',nlevel,' temp=',temp,' gam=',gamma_s,           &
             ' dzinc=',dzinc,' tb8=',tb8_k
      STOP
    END IF
    nlm1   = nlevel
    IF(nlm1 < 1) nlm1=1
    zht    = zbase
!
    DO j1=1,nlm1
      rl   = lathv+(273.15-temp)*dlvdt    ! Lv as func of Temp.
      arg  = rl*(temp-273.15)/273.15/temp/rv
      es   = eso*EXP(arg)                 ! satur. vapor press.
      arg  = -g*dzinc/rd/temp
      p    = press*EXP(arg)               ! hydrosta. press.decrease
!
!-----------------------------------------------------------------------
!
!  Calculate saturated adiabatic lapse rate
!
!-----------------------------------------------------------------------
!
      des   = es*rl/temp/temp/rv
      dtz   = -g*((1.0+0.621*es*rl/(press*rd*temp))/                    &
               (cp+0.621*rl*des/press))
      IF(ABS(dtz) < 1.0E-15) THEN
        PRINT*,' Zero dt/dz:',dtz,' g=',g,' es=',es,' press=',press
        PRINT*,'temp=',temp,' cp=',cp,' rl=',rl,' des=',des
        STOP
      END IF
      zht   = zht+dzinc  ! moist adiabatic ascent every 100m
      temp  = temp+dtz*dzinc
      IF(temp <= tb8_k) THEN   ! Cloud top is found
        ztop = zht + (tb8_k-temp)/dtz
        GO TO 150
      END IF
      press = p
    END DO  ! j=1,nlm1
    ztop = zht + (tb8_k-temp)/dtz
    150     CONTINUE
!
!-----------------------------------------------------------------------
!
!  Apply quality check to the cloud top height.
!
!-----------------------------------------------------------------------
!
    IF(ztop <= zs_1d(1)) THEN
      PRINT*,' Error, ztop for cloud top model is out of'               &
             ,' bound at grid pt(',i,j,').'
      PRINT*,' ztop=',ztop,' topo=',topo,                               &
             ' z_bottom=',zs_1d(1)
      PRINT*,'aborting...'
      STOP
    END IF
!
    IF(ztop > zs_1d(nz-1)) THEN
      PRINT*,' Warning, ztop for cloud top model is out of'             &
               ,' bound at grid pt(',i,j,').'
      PRINT*,' ztop=',ztop,' topo=',topo                                &
               ,' z_top=',zs_1d(nz-1)
      PRINT*,' ztop is set to the highest model level'
      cldtop_m_tb8 = zs_1d(nz-1)
      GO TO 900
    END IF
    cldtop_m_tb8 = ztop
    GO TO 900
!
    300     CONTINUE
!
!-----------------------------------------------------------------------
!
!  When Tb < -20 deg C, use the scheme that finds the cloud top
!  height by scanning the model temperature profile for a value
!  matching tb8.
!
!-----------------------------------------------------------------------
!
    tb8_loc=tb8_k
    IF(tb8_k <= tmin) THEN
      PRINT*,' Warning: cloud top colder than tmin in column'
      WRITE(6,*) '   i  j   tb8_k   t_top'
      WRITE(6,600) i,j,tb8_k,t_1d(nz-1)
      600       FORMAT(1X,2I3,2F8.2)
      PRINT*,' Cloud top temperature is set to tmin'
      tb8_loc=tmin
    END IF
!
    DO kl = nz-2,1,-1
!
!-----------------------------------------------------------------------
!
!  Find the lowest temp. crossing point in the model temperature
!  profile.
!
!-----------------------------------------------------------------------
!
      IF( (t_1d(kl)-tb8_loc) *  (t_1d(kl+1)-tb8_loc) <= 0.) THEN   ! Crossing Pt

        frac_k = (tb8_loc - t_1d(kl))                                   &
                    / (t_1d(kl+1) - t_1d(kl))
        arg = zs_1d(kl) + frac_k * (zs_1d(kl+1)-zs_1d(kl))

        IF(arg >= topo) THEN
          cldtop_m_tb8 = arg
        ELSE
          WRITE(6,*)' Warning: Cloud Top Below Ground - flagged'

          WRITE(6,601)
          601           FORMAT(1X,'  i  j kl  tb8_loc  t_abv  t_blw frac' &
                               ,'  hgt_m    topo  ctop_8 ')
          WRITE (6,602) i,j,kl,tb8_loc,t_1d(kl+1)                       &
                ,t_1d(kl),frac_k,arg,topo,cldtop_m_tb8
          602           FORMAT(1X,3I3,3F7.2,f5.2,3F8.0)
        END IF

      END IF   ! Crossing point found
    END DO ! kl

  ELSE ! No clouds according to SATELLITE (Band 8 - 11.2mm)
    l_tb8 = .false.

  END IF   ! tb8_k-t_gnd_k .lt. -thresh2
!
!-----------------------------------------------------------------------
!
!  Set variables
!
!-----------------------------------------------------------------------
!
  900   CONTINUE
  l_cloud_present = l_tb8
  cldtop_m = cldtop_m_tb8
  sat_cover = 1.0

  RETURN
END SUBROUTINE cloud_top





!
!##################################################################
!##################################################################
!######                                                      ######
!######             SUBROUTINE CORRECT_COVER                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE correct_cover(cover_in,cover_new_f                           & 1,2
           ,cldtop_old,cldtop_new,i,j,nx,ny,nz                          &
           ,zs_1d,t_1d,tb8_k,t_gnd_k                                    &
           ,istatus)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Find a thinner value for cloud cover consistent with the new
!  higher cloud top and the known brightness temperature.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Dan Birkenheuer?)
!  09/95.
!
!  MODIFICATION HISTORY:
!
!  04/29/96 (J. Zhang)
!  Modified for the ARPSDAS format.
!  Added full documentation.
!
!-----------------------------------------------------------------------
!
!  Variable Declaration
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!  INPUT:
  INTEGER :: i,j         ! grid point index
  INTEGER :: nx,ny,nz    ! grid size
  REAL :: t_1d(nz)       ! temp. analysis field
  REAL :: zs_1d(nz)      ! phycisal heights (m) for a column of scalar
                         ! grid points

  REAL :: tb8_k          ! satellite band 8 brightness temp.
  REAL :: t_gnd_k        ! ground skin temp.
  REAL :: cover_in       ! cloud cover before correction
  REAL :: cldtop_old     ! cloud top height before correction
!
!  OUTPUT:
  INTEGER :: istatus
  REAL :: cover_new_f    ! cloud cover after correction
  REAL :: cldtop_new     ! cloud top height after correction
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  REAL :: z_temp,temp_old,temp_new,frac
  REAL :: band8_cover    ! a function
  REAL :: cover_old,cover_new
  INTEGER :: iz_temp
  INTEGER :: iwrite
  DATA iwrite /0/
  SAVE iwrite
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Find Temperature of old cloud top
!
!-----------------------------------------------------------------------
!
  CALL hgt_to_zk (cldtop_old,z_temp,nz,zs_1d,istatus)
  IF(istatus /= 1)THEN
    PRINT*,' Old cldtop is out of domain range:', z_temp
    PRINT*,' Old cldtop=',cldtop_old,' zs range from',zs_1d(1)          &
            ,' to ',zs_1d(nz-1)
  END IF

  z_temp = MAX(1.,MIN(z_temp,FLOAT(nz-1)-.001))
  iz_temp = INT(z_temp)
  frac = z_temp - iz_temp
  temp_old = t_1d(iz_temp)    * (1. - frac)                             &
          +  t_1d(iz_temp+1)  * frac
!
!-----------------------------------------------------------------------
!
!  Find Temperature of new cloud top
!
!-----------------------------------------------------------------------
!
  CALL hgt_to_zk (cldtop_new,z_temp,nz,zs_1d,istatus)
  IF(istatus /= 1)THEN
    PRINT*,' New cldtop is out of domain range:', z_temp
    PRINT*,' New cldtop=',cldtop_new,' zs range from',zs_1d(1)          &
            ,' to ',zs_1d(nz-1)
!        return
  END IF

  z_temp = MAX(1.,MIN(z_temp,FLOAT(nz-1)-.001))
  iz_temp = INT(z_temp)
  frac = z_temp - iz_temp
  temp_new = t_1d(iz_temp) * (1. - frac)                                &
              + t_1d(iz_temp+1) * frac
!
!-----------------------------------------------------------------------
!
!  This one utilizes a linear approximation to
!  the sigma T**4 relationship
!
!-----------------------------------------------------------------------
!
  cover_old = MIN(cover_in,1.0)
  cover_new = cover_old * (tb8_k-t_gnd_k)/(temp_new-t_gnd_k)
!
!-----------------------------------------------------------------------
!
!  This one utilizes the sigma T**4 relationship
!
!-----------------------------------------------------------------------
!
  cover_new_f = band8_cover(tb8_k,t_gnd_k,temp_new)

  IF((j-1) == INT(j-1)/10*10) THEN
    iwrite = iwrite + 1
    IF(iwrite < 15) THEN
      WRITE(6,601,ERR=2)
      601       FORMAT(//1X,'**CORR_CVR**  i  j  tg_k   told   tnew '   &
                          ,'  ctold   ctnew   cvnew  cvrnewf')
      WRITE(6,1,ERR=2) i,j,t_gnd_k,temp_old,temp_new,cldtop_old         &
                           ,cldtop_new,cover_new,cover_new_f
      1         FORMAT(1X,12X,2I3,3F7.0,2F8.0,2F8.2)
    2       END IF
  END IF

  istatus=1

  RETURN
END SUBROUTINE correct_cover




!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION BAND8_COVER                    ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION band8_cover(tb8_k,t_gnd_k,t_cld)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Convert the radiance to the brightness temperature
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!  07/95
!
!  MODIFICATION HISTORY:
!  05/01/96  (Jian Zhang)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!  INPUT:
  REAL :: tb8_k       ! satellite band 8 brightness temp.
  REAL :: t_gnd_k     ! ground skin temp.
  REAL :: t_cld       ! estimated brightness temp. with SAO cld
!
!  OUTPUT:
  REAL :: band8_cover ! modified cld cover using sat. band 8 data
!
!  LOCAL:
  REAL :: r_sfc,r_sat,r_cld
  REAL :: temp_to_rad       ! a function
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  r_sfc = temp_to_rad(t_gnd_k)
  r_sat = temp_to_rad(tb8_k)
  r_cld = temp_to_rad(t_cld)

  band8_cover = (r_sat-r_sfc) / (r_cld-r_sfc)

  RETURN
  END FUNCTION band8_cover


!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION TEMP_TO_RAD                    ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION temp_to_rad(temp),1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Convert the radiance to the brightness temperature
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!  07/95
!
!  MODIFICATION HISTORY:
!  05/01/96  (Jian Zhang)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!  INPUT:
  REAL :: temp
!
!  OUTPUT:
  REAL :: temp_to_rad
!
!  LOCAL:
  INTEGER :: init
  DATA init /0/
  SAVE init

  INTEGER :: nsat
!
!  FUNCTIONS:  (in /vortex/lib/goeslib/invers.f)
  REAL :: vplanc       ! function that convert a brightness temperature
                     ! to the radiance for a specific band
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  INCLUDE 'adas.inc'

  IF(init == 0)THEN
    init = 1
    nsat = 3
    PRINT*,' Calling PLNKIV. calibration file:',calib_fname
    CALL plnkiv(nsat,calib_fname)
    PRINT*,' End calling PLNKIV'
  END IF

  temp_to_rad = temp
  temp_to_rad = vplanc(temp,8)

  RETURN
  END FUNCTION temp_to_rad


!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION RAD_TO_TEMP                    ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION rad_to_temp(rad),1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Convert the radiance to the brightness temperature
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!  07/95
!
!  MODIFICATION HISTORY:
!  05/01/96  (Jian Zhang)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!  INPUT:
  REAL :: rad
!
!  OUTPUT:
  REAL :: rad_to_temp
!
!  LOCAL:
  INTEGER :: init
  DATA init /0/
  SAVE init

  INTEGER :: nsat

!
!  FUNCTIONS:  (in /vortex/lib/goeslib/invers.f)
  REAL :: vbrite       ! function that convert a radiance to
                     ! a brightness temp for a specific band
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  INCLUDE 'adas.inc'

  IF(init == 0) THEN
    init = 1
    nsat = 3
    CALL plnkiv(nsat,calib_fname)
  END IF

  rad_to_temp = rad
  rad_to_temp = vbrite(rad,8)

  RETURN
  END FUNCTION rad_to_temp





!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE COMPARE_RAD                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE compare_rad (nx,ny,nz,clddiag,r_missing,zs,t_3d              & 1,6
           ,cldcv,t_sfc_k,t_gnd_k,tb8_k                                 &
           ,cvr_snow,nlyr)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  This routine compares the analyzed clouds to the 11.2mm radiation
!  and determines adjusted cloud fractions of the cloud layers to yield
!  a better fit.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (?)
!  09/95.
!
!  MODIFICATION HISTORY:
!
!  05/01/96 J. Zhang
!           Modified for the ARPSDAS format.
!           Added full documentation.
!  05/01/99 K. Brewster
!           Added check for missing IR data.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  INPUT:
  INTEGER :: nx,ny,nz     ! grid size
!
!  INPUT:
  INTEGER :: clddiag      ! print diagnostic fields?
  REAL :: r_missing       ! bad flag
  REAL :: t_3d(nx,ny,nz)  ! temp. analysis on the model grid
  REAL :: zs(nx,ny,nz)    ! physical heights (m) at model scalar pts.

  REAL :: cvr_snow(nx,ny) ! sfc snow coverage fraction
  REAL :: tb8_k(nx,ny)    ! satellite band8 brightness temp.
  REAL :: t_gnd_k(nx,ny)  ! ground skin temp.
  REAL :: t_sfc_k(nx,ny)  ! surface air temperature
!
!  OUTPUT:
  REAL :: cldcv(nx,ny,nz)   ! adjusted cloud cover analysis
!
!  LOCAL:
  REAL :: zs_1d(nz)             ! phy. heights (m) in a column of grid
  REAL :: t_1d(nz)              ! temp in a column of grid
  REAL :: cldcv_1d(nz)      ! cloud cover in one culomn of grid.
  REAL :: cldcv_1dref(nz)   ! cloud cover in one culomn of grid.

  REAL :: a(nz)                ! max. cldcvr in each cloud layer
  REAL :: f(nz)
  REAL :: a_new(nz)
  REAL :: f_new(nz)

  INTEGER :: ilyr(nz)       ! Dims needs to be changed to nz
  INTEGER :: ilyr_new(nz)   ! Dims needs to be changed to nz

  INTEGER :: nlyr,nlyr_new
  LOGICAL :: l_correct
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  REAL :: corr_thr,cover_step
  REAL :: delta_cover,delta_cover_ref
  INTEGER :: i_correct,iter,iter_max,n_iter
  REAL :: t_effective,tdiff,tdiff_ref
  INTEGER :: i,j,k,l,iwrite,imid,jmid,kmid,ibeg,iend
  INTEGER :: n_clear,n_clear_ns,n_clear_sn,n_cloudy,n_total             &
          ,n_correct
  REAL :: tdiff_sumsq,tdiff_cld_sumsq,tdiff_corr_sumsq                  &
      ,tdiff_corr_sumsq2,tdiff_corr_sumsq3,tdiff_corr_cld_sumsq
  REAL :: tb8_g_clr_sum,tb8_a_clr_sum,tb8_g_clr_sumsq                   &
      ,tb8_a_clr_sumsq,tb8_g_clr_ns_sum,tb8_a_clr_ns_sum                &
      ,tb8_g_clr_ns_sumsq,tb8_a_clr_ns_sumsq,tb8_g_clr_sn_sum           &
      ,tb8_a_clr_sn_sum,tb8_g_clr_sn_sumsq,tb8_a_clr_sn_sumsq

  REAL :: frac,frac_clouds
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  imid=nx/2
  jmid=ny/2
  ibeg=MAX( 1,imid-5)
  iend=MIN(nx,imid+5)
  kmid=nz/2
  IF(clddiag == 1) THEN
    WRITE(6,'(a)')' Comparing radiation'
    WRITE(6,611) nx,ny,nz
    611     FORMAT(1X,' nx=',i3,' ny=',i3,' nz=',i3)
!    612     FORMAT(1X,10F7.0)
    WRITE(6,'(a)')'grid point heights:'
    WRITE(6,613) (zs(i,jmid,kmid),i=ibeg,iend)
    613     FORMAT(10F8.0)
    WRITE(6,'(a)')'grid point temp:'
    WRITE(6,613) (t_3d(i,jmid,kmid),i=ibeg,iend)
    WRITE(6,'(a)')
!
    WRITE(6,'(a)') 'Surface temperatures: j=jmid, i=imid-5,+5'
    WRITE(6,614) (t_sfc_k(i,jmid),i=ibeg,iend)
    614     FORMAT(1X,10F7.0)
    WRITE(6,'(a)') 'Ground skin temperatures: j=jmid, i=imid-5,+5'
    WRITE(6,614) (t_gnd_k(i,jmid),i=ibeg,iend)
    WRITE(6,'(a)') 'Band 8 temperatures: j=jmid, i=imid-5,+5'
    WRITE(6,614) (tb8_k(i,jmid),i=ibeg,iend)
    WRITE(6,'(a)')
  END IF

  corr_thr = 4.             ! if diff of tb8-tb8e > 4K, correction
  cover_step = 0.05         ! add/minus 0.05 cldcvr each correction
  iter_max = 10             ! we'll do maximum of 10 iterns

  iwrite = 0                ! output message flag
  n_clear = 0               !
  n_clear_ns = 0
  n_clear_sn = 0
  n_cloudy = 0
  n_total = 0
  n_correct = 0
  n_iter = 0
  tdiff_sumsq = 0.          ! sum of sqr(tb8o-tb8e) for all pts
                             ! before correction
  tdiff_cld_sumsq = 0.      ! sum of sqare(tb8o-tb8e) for cloudy pts
                             ! before correction
  tdiff_corr_sumsq3 = 0.    ! sum of sqr(tb8o-tb8e) for cloudy pts
                            ! to be corrected, before the correction
  tdiff_corr_sumsq = 0.     ! sum of sqr(tb8o-tb8e) for all pts
                             ! after correction
  tdiff_corr_cld_sumsq = 0. ! sum of sqr(tb8o-tb8e) for all cldy pts
                             ! after correction
  tdiff_corr_sumsq2 = 0.    ! sum of sqr(tb8o-tb8e) for cloudy pts
                            ! to be corrected, after the correction
  tb8_g_clr_sum   = 0.      ! sum of tb8-tg for clear air pts
  tb8_a_clr_sum   = 0.      ! sum of tb8-tsfc for clear air pts
  tb8_g_clr_sumsq = 0.      ! sum of sqare(tb8-tg) for clear air pts
  tb8_a_clr_sumsq = 0.      ! sum of sqare(tb8-tsfc) for clear air pts
  tb8_g_clr_ns_sum   = 0.   ! for clear air with < 0.5 snow cover pts
  tb8_a_clr_ns_sum   = 0.
  tb8_g_clr_ns_sumsq = 0.
  tb8_a_clr_ns_sumsq = 0.
  tb8_g_clr_sn_sum   = 0.   ! for clear air with > 0.5 snow cover pts
  tb8_a_clr_sn_sum   = 0.
  tb8_g_clr_sn_sumsq = 0.
  tb8_a_clr_sn_sumsq = 0.

!
!-----------------------------------------------------------------------
!
!  Compare and correct (if necessary) the cloud cover analysis field
!  for each grid point.
!
!-----------------------------------------------------------------------
!
  DO j = 1,ny-1
    DO i = 1,nx-1
      IF(tb8_k(i,j) > 0.) THEN
!jz       do k = 1,nz-1
        DO k = 1,nz
          zs_1d(k) = zs(i,j,k)
          t_1d(k) = t_3d(i,j,k)
        END DO

        DO k = 1,nz
          cldcv_1dref(k) = cldcv(i,j,k)
        END DO

        CALL cvr_to_tb8e (nz,zs_1d,t_1d                                 &
                  ,cldcv_1dref,t_gnd_k(i,j)                             &
                  ,a,f,ilyr,t_effective,nlyr)
!
        tdiff = tb8_k(i,j)-t_effective ! Band 8 - calculated
        frac_clouds = 1.-f(nlyr)   ! column total cloud cover

        IF(nlyr == 1) THEN
          tb8_g_clr_sum   = tb8_g_clr_sum   + tdiff
          tb8_g_clr_sumsq = tb8_g_clr_sumsq + tdiff**2
          tb8_a_clr_sum   = tb8_a_clr_sum + tb8_k(i,j)-t_sfc_k(i,j)
          tb8_a_clr_sumsq=tb8_a_clr_sumsq+(tb8_k(i,j)-t_sfc_k(i,j))**2
          n_clear = n_clear + 1

          IF(cvr_snow(i,j) /= r_missing) THEN
            IF(cvr_snow(i,j) < 0.5) THEN
              tb8_g_clr_ns_sum   = tb8_g_clr_ns_sum   + tdiff
              tb8_g_clr_ns_sumsq = tb8_g_clr_ns_sumsq + tdiff**2
              tb8_a_clr_ns_sum = tb8_a_clr_ns_sum                       &
                                  + tb8_k(i,j)-t_sfc_k(i,j)
              tb8_a_clr_ns_sumsq = tb8_a_clr_ns_sumsq                   &
                                  + (tb8_k(i,j) - t_sfc_k(i,j))**2
              n_clear_ns = n_clear_ns + 1
            ELSE
              tb8_g_clr_sn_sum   = tb8_g_clr_sn_sum   + tdiff
              tb8_g_clr_sn_sumsq = tb8_g_clr_sn_sumsq + tdiff**2
              tb8_a_clr_sn_sum = tb8_a_clr_sn_sum                       &
                                  + tb8_k(i,j)-t_sfc_k(i,j)
              tb8_a_clr_sn_sumsq = tb8_a_clr_sn_sumsq                   &
                                  +(tb8_k(i,j)-t_sfc_k(i,j))**2
              n_clear_sn = n_clear_sn + 1
            END IF
          END IF
        ELSE
          n_cloudy = n_cloudy + 1
          tdiff_cld_sumsq = tdiff_cld_sumsq + tdiff**2
        END IF

        n_total = n_total + 1
        tdiff_sumsq = tdiff_sumsq + tdiff**2
!
!-----------------------------------------------------------------------
!
!    Apply corrections?
!
!    Only when the following four conditions are satisfied:
!    1. at least one cloud layer
!    2. the difference between the observed and estimated tb8
!       is larger than the threshold (4K)
!    3. tb8_obs - t_gnd < -8K
!    4. column total cloud cover is larger than 40%
!
!-----------------------------------------------------------------------
!
        IF (nlyr >= 2 .AND. ABS(tdiff) > corr_thr                       &
              .AND. t_gnd_k(i,j) - tb8_k(i,j) > 8.                      &
              .AND. frac_clouds > 0.4 )  THEN
          l_correct = .true.
          n_correct = n_correct + 1
          tdiff_corr_sumsq3 = tdiff_corr_sumsq3 + tdiff**2
        ELSE
          l_correct = .false.
        END IF
!
!-----------------------------------------------------------------------
!
!    This corrective section is now turned on
!
!-----------------------------------------------------------------------
!
        iter = 0
        delta_cover = 0.

        905       IF (l_correct) THEN
          iter = iter + 1
          tdiff_ref = tdiff             ! save initial difference
          delta_cover_ref = delta_cover ! save initial correction

          IF(tdiff < 0.) THEN
            delta_cover = delta_cover + cover_step
          ELSE
            delta_cover = delta_cover - cover_step
          END IF

          CALL apply_correction (cldcv_1dref,cldcv_1d                   &
                  ,nz,ilyr,f,delta_cover)

          CALL cvr_to_tb8e (nz,zs_1d,t_1d                               &
                  ,cldcv_1d,t_gnd_k(i,j)                                &
                  ,a_new,f_new,ilyr_new,t_effective,nlyr_new)
!
          tdiff = tb8_k(i,j)-t_effective
          frac_clouds = 1.-f_new(nlyr_new)
!
!-----------------------------------------------------------------------
!
!  Continue to apply corrections?
!
!-----------------------------------------------------------------------
!
          IF(nlyr_new >= 2 .AND. frac_clouds > 0.4) THEN
            i_correct = 1
          ELSE
            i_correct = 0
          END IF

          IF(MOD(iwrite,50) == 0) THEN
            iwrite = iwrite + 1
            WRITE(6,*) 'iter=',iter ,' tdiff_ref=',tdiff_ref            &
                       ,'tg=',t_gnd_k(i,j)
            WRITE(6,*) (a(l),l=nlyr-1,1,-1)
            WRITE(6,641)
            641           FORMAT(1X,' i  j  ncn  tb8o   tb8e  tdiff  cldcvm i_c' &
                                ,6(' cvrln'))
            WRITE(6,642,ERR=912)i,j,nlyr_new-1,tb8_k(i,j),t_effective   &
                 ,tdiff,frac_clouds,i_correct                           &
                 ,(a_new(l),l=nlyr_new-1,1,-1)
            642           FORMAT(1X,2I3,i4,2F7.0,f7.1,f7.3,i2,2X,10F6.2)
          912         END IF

          IF(i_correct == 1 .AND. iter < iter_max                       &
              .AND. ABS(tdiff) < ABS(tdiff_ref)                         &
              .AND. tdiff * tdiff_ref > 0.) GO TO 905
! Loop back & increment cover
          n_iter = n_iter + iter
!
!-----------------------------------------------------------------------
!
!      Final iteration
!
!-----------------------------------------------------------------------
!
          IF(.false. .OR. tdiff*tdiff_ref >= 0.) THEN
                                      ! Select best of last two increments
            IF(ABS(tdiff) >= ABS(tdiff_ref)) THEN
              delta_cover = delta_cover_ref
              CALL apply_correction (cldcv_1dref,cldcv_1d               &
                    ,nz,ilyr,f,delta_cover)
              tdiff = tdiff_ref
            END IF
          ELSE           ! Do one Newton iteration
            frac = tdiff / (tdiff - tdiff_ref)
            delta_cover=delta_cover+frac*(delta_cover_ref-delta_cover)

            CALL apply_correction (cldcv_1dref,cldcv_1d                 &
                  ,nz,ilyr,f,delta_cover)

            CALL cvr_to_tb8e (nz,zs_1d,t_1d                             &
                  ,cldcv_1d,t_gnd_k(i,j)                                &
                  ,a_new,f_new,ilyr_new,t_effective,nlyr_new)

            tdiff = tb8_k(i,j)-t_effective
            n_iter = n_iter + 1
          END IF
!
!-----------------------------------------------------------------------
!
!  Put the corrected column cloud coverback to 3D cloud cover field
!
!-----------------------------------------------------------------------
!
          DO k=1,nz
            cldcv(i,j,k) = cldcv_1d(k)
          END DO

        END IF ! Corrections were made

        tdiff_corr_sumsq     = tdiff_corr_sumsq + tdiff**2
        IF(nlyr > 1)      & ! cloudy (at least before adjustment)
        tdiff_corr_cld_sumsq = tdiff_corr_cld_sumsq + tdiff**2

        IF(l_correct) THEN
          tdiff_corr_sumsq2 = tdiff_corr_sumsq2 + tdiff**2
        END IF

      END IF

    END DO ! I
  END DO ! J
!
!-----------------------------------------------------------------------
!
!  Write out statistics on consistency of Band 8 data and cloud cover
!
!-----------------------------------------------------------------------
!
  IF(n_clear > 0) THEN
    WRITE(6,951,ERR=960) n_clear, tb8_g_clr_sum/FLOAT(n_clear)          &
                   ,SQRT(tb8_g_clr_sumsq/FLOAT(n_clear))
    951     FORMAT(/' Mean/RMS band 8/gnd temp residual in clear skies =' &
                   ,i5,2F9.3)
    960     WRITE(6,961,ERR=970) n_clear, tb8_a_clr_sum/FLOAT(n_clear)  &
                           ,SQRT(tb8_a_clr_sumsq/FLOAT(n_clear))
    961     FORMAT(' Mean/RMS band 8/air temp residual in clear skies =' &
                   ,i5,2F9.3)
  970   END IF

  IF(n_clear_ns > 0) THEN
    WRITE(6,956,ERR=965) n_clear_ns                                     &
                         ,tb8_g_clr_ns_sum/FLOAT(n_clear_ns)            &
                       ,SQRT(tb8_g_clr_ns_sumsq/FLOAT(n_clear_ns))
    956     FORMAT(/' Mean/RMS band 8/gnd temp resid in clear/nsnow skies =' &
                    ,i5,2F9.3)
    965     WRITE(6,966,ERR=975) n_clear_ns                             &
                                 ,tb8_a_clr_ns_sum/FLOAT(n_clear_ns)    &
                              ,SQRT(tb8_a_clr_ns_sumsq/FLOAT(n_clear_ns))
    966     FORMAT(' Mean/RMS band 8/air temp resid in clear/nsnow skies =' &
                    ,i5,2F9.3)
  975   END IF

  IF(n_clear_sn > 0) THEN
    WRITE(6,976,ERR=985) n_clear_sn                                     &
                         ,tb8_g_clr_sn_sum/FLOAT(n_clear_sn)            &
                     ,SQRT(tb8_g_clr_sn_sumsq/FLOAT(n_clear_sn))
    976     FORMAT(/' Mean/RMS band 8/gnd temp resid in clear/snow skies =' &
                    ,i5,2F9.3)
    985     WRITE(6,986,ERR=995) n_clear_sn                             &
                                 ,tb8_a_clr_sn_sum/FLOAT(n_clear_sn)    &
                              ,SQRT(tb8_a_clr_sn_sumsq/FLOAT(n_clear_sn))
    986     FORMAT(' Mean/RMS band 8/air temp resid in clear/snow skies =' &
                    ,i5,2F9.3)
  995   END IF

  IF(n_total > 0) THEN
    WRITE(6,971,ERR=980)n_total,SQRT(tdiff_sumsq/FLOAT(n_total))
    971     FORMAT(/' RMS band 8/teff residual (bfr corr) in all skies =' &
                    ,i5, f9.3)
    980     WRITE(6,981,ERR=990) n_total,                               &
                                 SQRT(tdiff_corr_sumsq/FLOAT(n_total))
    981     FORMAT(' RMS band 8/teff residual (aft corr) in all skies =' &
                    ,i5, f9.3)
  990   END IF

  IF(n_cloudy > 0) THEN
    WRITE(6,991,ERR=1000) n_cloudy,                                     &
                          SQRT(tdiff_cld_sumsq/FLOAT(n_cloudy))
    991     FORMAT(' RMS band 8/teff residual (bfr corr) in cldy skies =' &
                    ,i5, f9.3)
    1000    WRITE(6,1001,ERR=1010) n_cloudy,                            &
                                  SQRT(tdiff_corr_cld_sumsq/FLOAT(n_cloudy))
    1001    FORMAT(' RMS band 8/teff residual (aft corr) in cldy skies =' &
                    ,i5, f9.3)
  1010  END IF

  IF(n_correct > 0) THEN
    WRITE(6,1011,ERR=1020)                                              &
            n_correct,SQRT(tdiff_corr_sumsq3/FLOAT(n_correct))
    1011    FORMAT(/' RMS band 8/teff resid (bfr corr - corrected pts) =' &
                    ,i5, f9.3)
    1020    WRITE(6,1021,ERR=1030)                                      &
                    n_correct,SQRT(tdiff_corr_sumsq2/FLOAT(n_correct))
    1021    FORMAT(' RMS  band 8/teff resid (aft corr - corrected pts) =' &
                    ,i5, f9.3)
    WRITE(6,*)' Total/Average # of iterations = ',n_iter,               &
               FLOAT(n_iter)/FLOAT(n_correct)
    WRITE(6,*)
  1030  END IF

  RETURN
END SUBROUTINE compare_rad





!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE CVR_TO_TB8E                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cvr_to_tb8e (nz,zs_1d,t_1d                                   & 3
           ,cvr,t_gnd_k                                                 &
           ,a,f,ilyr,t_effective,nlyr)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  To calculate the estimated band8 brightness temperature using
!  the cloud cover analysis.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  (Jian Zhang)
!  04/1996   Based on the LAPS cloud analysis code of 09/95
!
!  MODIFICATION HISTORY:
!
!  04/26/96  J. Zhang
!            Modified for the ARPSDAS format. Added full
!            documentation.
!  07/27/96  J. Zhang
!            Modified code so that the cloud layer with top at the
!            upper boundary of the cloud grid is counted.
!  04/11/97  J. Zhang
!            Include adascld24.inc for ncloud
!  08/06/97  J. Zhang
!            Change adascld24.inc to adascld25.inc
!  09/11/97  J. Zhang
!            Change adascld25.inc to adascld26.inc
!  04/20/98  J. Zhang
!            Based on the arps4.3.3 version. Abandoned cloud grid,
!            using the model grid.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  INPUT
  INTEGER :: nz
!
!  INPUT:
  REAL :: cvr(nz)     ! Cloud cover from analysis

  REAL :: zs_1d(nz)
  REAL :: t_1d(nz)
  REAL :: t_gnd_k
!
!  OUTPUT:
  REAL :: t_effective    ! effective brightness temp. estimated
                         ! from cloud cover analysis
  INTEGER :: nlyr
  INTEGER :: ilyr(nz)    ! Layer index for each cloud lvl (needs nz)

  REAL :: a(nz)          ! Cloud fractions of layers
  REAL :: f(nz)          ! Apparnt "x-sectn" of cldlyrs seen from above
!
!  LOCAL:
  INTEGER :: ik(nz)      ! Height level representative of cloud layers
  REAL :: temp_lyr(nz)   ! temp. at "ik" lavels
!
!  FUNCTIONS:
  REAL :: temp_to_rad
  REAL :: rad_to_temp
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: k,n
  REAL :: rsum,sumf
  REAL :: cvr_thresh
  PARAMETER (cvr_thresh = 0.1)
!
!-----------------------------------------------------------------------
!
!  Convert from cld cvr to discreet cloud layer indices (cvr to a)
!
!  Find the maximum cloud cover "a" in each cloud deck, and
!  the level index "ik" associated with the maximum cloud cover.
!
!-----------------------------------------------------------------------
!
  nlyr = 0

  IF(cvr(nz) >= cvr_thresh) THEN  ! cldtop at the upper boundary
    nlyr = nlyr + 1
    a(nlyr) = cvr(nz)
    ik(nlyr) = nz
    ilyr(nz) = nlyr
  ELSE
    ilyr(nz)=0
  END IF

  DO k = nz-1,1,-1
    IF(cvr(k) >= cvr_thresh .AND. cvr(k+1) < cvr_thresh) THEN
      nlyr = nlyr + 1
      a(nlyr) = cvr(k)
      ik(nlyr) = k
    ELSE
      IF(nlyr >= 1) THEN
        IF(cvr(k) > a(nlyr)) THEN
          a(nlyr) = cvr(k)         ! Max cldcvr within a layer
          ik(nlyr) = k             ! the lvl index for the max. cldcv
        END IF
      END IF
    END IF

    IF(cvr(k) >= cvr_thresh) THEN      ! Still within layer
      ilyr(k) = nlyr
    ELSE                                 ! Below layer
      ilyr(k) = 0
    END IF

  END DO ! k
!
!-----------------------------------------------------------------------
!
!  Get temperatures of the maximum cldcvr level in each cld deck.
!
!-----------------------------------------------------------------------
!
  DO n = 1,nlyr
    k = ik(n)
    temp_lyr(n) = t_1d(k)
  END DO ! n

!
!-----------------------------------------------------------------------
!
!  Add a layer for the ground
!
!-----------------------------------------------------------------------
!
  nlyr = nlyr + 1
  a(nlyr) = 1.0
  temp_lyr(nlyr) = t_gnd_k
!
!-----------------------------------------------------------------------
!
!  Convert cloud layer fractions to "cross-section" seen from
!  satellite.  This solves for the "f" array given the "a" array
!
!-----------------------------------------------------------------------
!
  a(1) = MIN(  f(1) = a(1)
  sumf = f(1)
  IF(nlyr >= 2) THEN
    DO n = 2,nlyr
      a(n) = MIN(      f(n) = a(n) * (1.0 - sumf) ! fraction of radiation reaches
                                 ! top of atm from each cld layer
      sumf = sumf + f(n)         ! fraction of radiation blked
                                 ! /attened by all cld lyrs above
    END DO ! n
  END IF ! nlyr
!
!-----------------------------------------------------------------------
!
!  Calculate total radiance from all cloud layers + ground
!
!-----------------------------------------------------------------------
!
  rsum = 0
  DO n = 1,nlyr
    rsum = rsum + temp_to_rad(temp_lyr(n)) * f(n)
  END DO ! n
!
!-----------------------------------------------------------------------
!
!  Convert to effective temperature and compare to
!  the observed brightness temp
!
!-----------------------------------------------------------------------
!
  t_effective = rad_to_temp(rsum)

  RETURN
END SUBROUTINE cvr_to_tb8e


!
!##################################################################
!##################################################################
!######                                                      ######
!######            SUBROUTINE APPLY_CORRECTION               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE apply_correction (cldcv_1dref,cldcv_1d                       & 3
           ,nz,ilyr,f,delcv)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  To apply the correction to the cloud cover field.
!  The amount for correction is an input parameter.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  05/02/96
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  INPUT:
  INTEGER :: nz        ! number of cloud analysis levels
  REAL :: cldcv_1dref(nz) ! the column cldcvr to be corr.
  INTEGER :: ilyr(nz)  ! cld deck(layer) indices for ea.cldy level
  REAL :: f(nz)           ! fracn of rad. blocked by each cloud deck
  REAL :: delcv            ! the cld cvr increment added to the input
!
!  OUTPUT:
  REAL :: cldcv_1d(nz)     ! corrected cloud cover
!
!  LOCAL:
  INTEGER :: k
!
!-----------------------------------------------------------------------
!
!  Apply correction to 1D cloud cover field
!
!-----------------------------------------------------------------------
!
  DO k = 1,nz
!jz     if( ilyr(k).gt.0 .and. f(ilyr(k)).gt.0.0) then
    IF(ilyr(k) > 0) THEN
      IF(        cldcv_1d(k) = MIN(cldcv_1dref(k)+delcv, 1.0)
        cldcv_1d(k) = MAX(cldcv_1d(k), 0.0)
      ELSE
        cldcv_1d(k) = MAX(MIN(cldcv_1dref(k), 1.0), 0.0)
      END IF
    ELSE
      cldcv_1d(k) = MAX(MIN(cldcv_1dref(k), 1.0), 0.0)
    END IF
  END DO

  RETURN
END SUBROUTINE apply_correction