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


SUBROUTINE barnes_r5 (nx,ny,nz,TO,t,wt_p,cf_modelfg,                    & 1
           l_stn_clouds,default_clear_cover,cld_snd,wt_snd,             &
           spcng,i_snd,j_snd,n_cld_snd,                                 &
           istatus)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Horizontally interpolate the SAO cloud soundings to ADAS grid.
!  The interpolation scheme is Barnes-type, the weight given to a
!  data point is:
!
!  wt(m) = [(100.0*bias)/(1.533*dgrd(m)**2 + 1.0)]**(5.0/2.0)
!  wt(m) ~ 34367.237 /(dgrd(m)**5)
!
!  where
!
!    m  is the index of a data point
!    bias=1.0        ! it can be a function of data location
!    dgrd(m) is the distance from the mth data point to
!          the grid point, in terms of the number of grid spacings.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  (Jian Zhang)
!  07/95   Based on the LAPS cloud analysis code
!
!  MODIFICATION HISTORY:
!  02/06/96  J. Zhang
!            Modified for ADAS. Added documents.
!  05/10/96  J. Zhang
!            Added the interpolation for the station grid column
!            (skipped during the first step).
!  03/14/97  J. Zhang
!            Clean up the code and implemented for the offcial
!            arps4.2.4 version
!  07/30/97  J. Zhang
!            Added minimum weight threshold for Barnes interpolation.
!            This can avoid the interpolations at the grid points
!            where there is no station nearby.
!  08/06/97  J. Zhang
!            Change adascld24.inc to adascld25.inc.
!            Set weight_modelfg=1.0 for the cases when there is
!            no SAO cloud soundings.
!  09/09/97  J. Zhang
!            Using the different influence cutoff radius for
!            cloud soundingd with different cloud coverage
!  09/10/97  J. Zhang
!            Awoke the procedure which uses background RH field to
!            derive cloud cover field when there is no SAO cloud
!            sounding.
!            Change adascld25.inc to adascld26.inc.
!  11/18/97  J. Zhang
!            Moved "max_obs" to adascld26.inc.
!  11/19/97  J. Zhang
!            Added test for the case when ncnt > max_obs.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!
!  OUTPUT:
!
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'adas.inc'
!
!   INPUT from adas.inc:
!   real r_missing      ! missing value flag
!   real range_cld      ! cutoff radius for interpolation
!   integer max_cld_snd ! max. possible # of cloudy soundings
!   integer max_obs     ! max. possible # of cldy data points in 3D
                        ! domain (should be about max_cld_snd*nz)
!   integer i_perimeter ! an extension of the ADAS domain for
                        ! searching SAOs.
!
! INPUT:
!
  INTEGER :: nx,ny,nz           ! grid size
  INTEGER :: n_cld_snd          ! actual # of cld snd

  REAL :: TO(nx,ny,nz)          ! 1st guess of gridded cloud cover
  REAL :: cf_modelfg(nx,ny,nz)  ! 1st guess for cloud cover field
  REAL :: wt_p(nx,ny,nz)        ! weights for first guess field

  LOGICAL :: l_stn_clouds       ! using SAO stns'cld sndings?
  REAL :: default_clear_cover   ! defaule cld value for clear sky

  REAL*4 cld_snd(max_cld_snd,nz) ! obs. cloud cover sounding
  REAL*4 wt_snd(max_cld_snd,nz)  ! weights for obs. cld sounding

  INTEGER*4 i_snd(max_cld_snd)   ! i-loc of cloud sounding stn
  INTEGER*4 j_snd(max_cld_snd)   ! j-loc of cloud sounding stn
!
! OUTPUT:
! 
  REAL :: t(nx,ny,nz)       ! the gridded analysis of cld cvr
  INTEGER :: istatus
!
! LOCAL:
!
  INTEGER :: nx_lut,ny_lut      ! dims for look up tbl of intp wgt.
  INTEGER :: ix_low,ix_high,iy_low,iy_high

  INTEGER :: n_fnorm
!
  INTEGER*4 iskip        ! # of grid pts skipped when performing
                         ! Barnes intp. (for time-saving)
  INTEGER*4 lowi_lut(nx) ! index for skipped grid points.
  INTEGER*4 lowj_lut(ny) ! index for skipped grid points.
!
  REAL, allocatable :: fnorm(:)    ! normalized weights
!
  INTEGER :: iob(max_obs)   ! i-loc of each cld data pt. in ADAS grid
  INTEGER :: job(max_obs)   ! j-loc of each cld data pt. in ADAS grid
  INTEGER :: kob(max_obs)   ! k-lvl of each cld data pt. in ADAS grid
  INTEGER :: nob(max_obs)   ! sequential index of each cldy data pt.

  REAL, allocatable :: iiilut(:,:)
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nlast (nz)        ! the index for the cloud sounding data
                             ! point at each level.
  LOGICAL :: l_analyze (nz)
!
  REAL :: weight_modelfg       ! the weight given to the 1st guess
                             ! cloud cover field
  REAL :: rden_ratio,radm2,rfac
  PARAMETER ( rfac=10.0, radm2=1.533/rfac**2)
!   parameter ( radm2=1.533)
!
  REAL :: exp_dist_wt          ! power of distance in weight
  PARAMETER( exp_dist_wt = 5.0)

  REAL :: spcng                ! grid spacing in m

  REAL :: rr,rr_max
  INTEGER :: iii,iii_max

  INTEGER :: iiizero,ncnt
  REAL :: bias_iii             ! a bias may be a func of locatns

  LOGICAL :: cldprt
  INTEGER :: i,j,k,n,nobs,nstart,nstop,ii,jj,nn,n1,jmid,ibeg
  REAL :: weight,sum,sumwt

  INTEGER :: il,ih,jl,jh,km1
  REAL :: z1,z2,z3,z4,fraci,fracj
  INTEGER :: astat
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  nx_lut = nx + i_perimeter - 1
  ny_lut = ny + i_perimeter - 1
  ix_low = 1  - i_perimeter
  ix_high= nx + i_perimeter
  iy_low = 1  - i_perimeter
  iy_high= ny + i_perimeter
  n_fnorm = 1.6 * ((nx_lut*nx_lut)+(ny_lut*ny_lut))

  allocate (fnorm(n_fnorm),stat=astat)
  IF (astat /= 0) THEN
    WRITE (6,*) "BARNES_R5: ERROR allocating fnorm, exiting"
    STOP
  END IF
  allocate (iiilut(-nx_lut:nx_lut,-ny_lut:ny_lut),stat=astat)
  IF (astat /= 0) THEN
    WRITE (6,*) "BARNES_R5: ERROR allocating fnorm, exiting"
    STOP
  END IF

  cldprt=(clddiag == 1)
  jmid=ny/2
  ibeg=MAX(1,nx-10)
  WRITE(6,1234)
  1234  FORMAT(1X,'Barnes_r5 called')

  istatus = 0
  IF(cldprt) THEN
    PRINT*, ' # of cloud soundings = ',n_cld_snd                        &
         ,'  wfg=',weight_modelfg
!
    WRITE(6,'(/a)') ' ==== BARNES cloud cover first guess===='
    WRITE(6,632) (i,i=nx-10,nx)
    632     FORMAT(/1X,3X,11(1X,i4,1X))
    WRITE(6,633) (k,(cf_modelfg(i,jmid,k),i=ibeg,nx)                    &
                ,k=nz,1,-1)
    633     FORMAT(1X,i3,11F6.1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Obtain and/or iterate for value of iskip
!
!-----------------------------------------------------------------------
!
  iskip = nint(20000. / spcng)           ! # of grid pts skipped
                                          ! during calcul'n
  iskip = MAX(iskip,1)

  100   rden_ratio = ((FLOAT(nx)-1.)/FLOAT(iskip))

  IF( ABS(rden_ratio - nint(rden_ratio)) > 0.001) THEN
    WRITE(6,*)' Bad value of iskip'
    IF(iskip > 1)THEN
      iskip = iskip - 1
      GO TO 100
    ELSE IF( iskip == 1) THEN
      WRITE(6,*)' Code error - stop'
      STOP
    END IF
  ELSE
    WRITE(6,*)' Good value of iskip = ',iskip
  END IF
!
!-----------------------------------------------------------------------
!
!  Create an array (lookup table) of weights as function of the
!  distance iii (in the unit of grid spacings)
!
!-----------------------------------------------------------------------
!
  IF( l_stn_clouds) PRINT*, ' # of cloud soundings = ',n_cld_snd

  iiizero = 0

  DO iii = 1,n_fnorm
    bias_iii = 1E0    ! It can be func of location iii later

    fnorm(iii) = (100.*bias_iii/FLOAT(iii)) ** (exp_dist_wt/2.0)
    IF( fnorm(iii) == 0. .AND. iiizero == 0) THEN
      iiizero = iii
      WRITE(6,*)' WARNING: fnorm array reached zero, iii=',iiizero
    END IF
  END DO  ! III

  WRITE(6,*)' Min possible weight = ',fnorm(n_fnorm)
  WRITE(6,*)' Max possible weight = ',fnorm(1)

  IF (n_cld_snd <= 0) THEN
    weight_modelfg = 1.0
  ELSE
    weight_modelfg = 0.0
  END IF
!
!-----------------------------------------------------------------------
!
!  Count the number of cloudy data points on each cloud height
!  level, and mark each cloudy data point with iob- and job-index.
!
!-----------------------------------------------------------------------
!
  ncnt=0        ! # of cloudy data points on each cloud grid level

  DO k = 1, nz

    IF (.NOT. l_stn_clouds) THEN   ! using gridded non-model 1st
                                   ! guess field
      nlast(k) = ncnt
      DO j=1,ny
        DO i=1,nx
          IF( TO(i,j,k) == r_missing) GO TO 223
          ncnt=ncnt+1
          iob(ncnt)=i
          job(ncnt)=j
          kob(ncnt)=k
          nlast(k) = ncnt
          223         CONTINUE
        END DO ! i
      END DO ! j

    ELSE ! l_stn_clouds = .TRUE.,  use cloud soundings
      nlast(k) = ncnt

      IF (n_cld_snd >= 1) THEN
        DO n = 1,n_cld_snd
          IF(cld_snd(n,k) < default_clear_cover) GO TO 233
!
!-----------------------------------------------------------------------
!
!  Test if out of the bounds of the ADAS_plus domain
!
!-----------------------------------------------------------------------
!
          IF(i_snd(n) < ix_low .OR. i_snd(n) > ix_high) GO TO 233
          IF(j_snd(n) < iy_low .OR. i_snd(n) > iy_high) GO TO 233

          ncnt=ncnt+1
          IF (ncnt > max_obs) THEN
            PRINT*,'# of cldy grid pts exceeds the max. # allowed'
            PRINT*,'# of cldy grid pts:',ncnt,' # allowed:',max_obs
            PRINT*,'Plz increase MAX_OBS in adas.inc.'
            PRINT*,' Aborting...'
            STOP
          END IF
          iob(ncnt)=i_snd(n)
          job(ncnt)=j_snd(n)
          kob(ncnt)=k
          nob(ncnt)=n
          nlast(k) = ncnt
          233           CONTINUE
        END DO ! n
      END IF

    END IF ! l_stn_clouds

    IF( k > 1) THEN
      km1 = k - 1
      l_analyze(k) = .false.

      IF(.NOT. l_stn_clouds) THEN
        DO j=1,ny
          DO i=1,nx
            IF( TO(i,j,k) /= TO(i,j,km1)) THEN
              l_analyze(k) = .true.
              GO TO 250
            END IF
          END DO ! i
        END DO ! j

      ELSE
        DO n=1,n_cld_snd
          IF( cld_snd(n,k) /= cld_snd(n,km1)) THEN
            l_analyze(k) = .true.
            goto 250
          END IF
        END DO ! n

      END IF ! l_stn_clouds

      250       CONTINUE

    ELSE ! k .eq. 1
      l_analyze(k) = .true.
    END IF

  END DO ! K

  IF(ncnt == 0) THEN
    PRINT*, 'No SAO cloud sounding data for Barnes'
    PRINT*, 'Using the background cloud cover field'

    DO k=2,nz-2
      DO j=1,ny
        DO i=1,nx
          t(i,j,k) = cf_modelfg(i,j,k)
        END DO
      END DO
    END DO
    istatus = 1
    RETURN
  ELSE
    WRITE(6,*)' Ncnt = ',ncnt
  END IF
!
!-----------------------------------------------------------------------
!
!  Create a lookup table for fnorm(iii): the weights as a function
!  of i- and j- distances.
!
!-----------------------------------------------------------------------
!
  rr_max = (nx_lut-1)**2 + (ny_lut-1)**2
  iii_max = radm2*100.*rr_max + 1.
  PRINT*,' radm2*100,iii_max,n_fnorm: ',radm2*100.,iii_max,n_fnorm

  IF( iii_max > n_fnorm) THEN
    PRINT*,' iii_max is too large, increase n_fnorm',iii_max            &
          ,n_fnorm
    istatus = 0
    RETURN
  END IF
!
  DO i = -nx_lut,nx_lut
    DO j = -ny_lut,ny_lut
      rr=i*i+j*j
      iii=radm2*100.*rr+1.
      IF( iii > n_fnorm) iii=n_fnorm
      iiilut(i,j) = fnorm(iii)            ! lookup table for weights
                                          ! fnorm as a function of
                                          ! distance iii.
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Loop through each cloud height level.
!
!-----------------------------------------------------------------------
!
  DO k = 1, nz

    IF(k == 1)THEN
      nstart = 1
    ELSE
      nstart = nlast(k-1) + 1
    END IF

    nstop = nlast(k)
    nobs = nstop-nstart+1       ! # of obs. in one level

    IF (nobs <= 0) THEN  ! No Obs; Set level to model 1st guess
!
      IF(cldprt) WRITE(6,42) k,nstart,nstop,nobs
      42        FORMAT(' lvl,nstart,nstop,nobs=',4I5,                   &
                              ' No Obs; Set level to model fg')
      DO j=1,ny
        DO i=1,nx
          t(i,j,k) = cf_modelfg(i,j,k)
        END DO ! i
      END DO ! j

    ELSE IF( weight_modelfg > 1.0E-15 .OR. l_analyze(k) ) THEN

      IF(cldprt) WRITE(6,50) k,nstart,nstop,nobs
      50        FORMAT(' lvl,nstart,nstop,nobs=',4I5)
!
!-----------------------------------------------------------------------
!
!    Analyze every iskip-th  grid point
!
!-----------------------------------------------------------------------
!
      DO j=1,ny,iskip
        DO i=1,nx,iskip
          sum=0.
          sumwt=0.
          IF(.NOT. l_stn_clouds)THEN
            DO n=nstart,nstop
              ii=iob(n)
              jj=job(n)
              iii = (i-ii)*(i-ii) + (j-jj)*(j-jj)
              rr = FLOAT (iii)
              rr = spcng*SQRT (rr)
              IF (rr <= range_cld) THEN
!
!-----------------------------------------------------------------------
!
!c              gridded obs are being weighted
!
!-----------------------------------------------------------------------
!
                weight = iiilut(i-ii,j-jj) * wt_p(ii,jj,k)
                sum=weight*TO(ii,jj,k)+sum
                sumwt=sumwt+weight
              END IF
            END DO  ! n
          ELSE
            DO n=nstart,nstop
              ii=iob(n)
              jj=job(n)
              nn=nob(n)
              iii = (i-ii)*(i-ii) + (j-jj)*(j-jj)
              rr = FLOAT (iii)
              rr = spcng*SQRT (rr)
              IF (rr <= range_cld) THEN
!
!-----------------------------------------------------------------------
!
!c              Station obs are being weighted
!
!-----------------------------------------------------------------------
!
                weight = iiilut(i-ii,j-jj) * wt_snd(nn,k)
                sum=weight*cld_snd(nn,k)+sum
                sumwt=sumwt+weight
              END IF
            END DO ! n
          END IF

          sum   = sum   + weight_modelfg * cf_modelfg(i,j,k)
          sumwt = sumwt + weight_modelfg

          IF (ABS(sumwt) <= 1.0E-15) THEN
!9/3/97        t(i,j,k) = r_missing
!              istatus = 0
!              print*,'(00): ',i,j,k,' tot sum:',sum,'  tot sumwt:'
!     :                ,sumwt,'  t:',t(i,j,k)
            t(i,j,k) = 0.01
          ELSE
            t(i,j,k)=sum/sumwt
          END IF

        END DO ! i
      END DO ! j
!
!-----------------------------------------------------------------------
!
!      Bilinearly interpolate to fill in rest of domain
!
!-----------------------------------------------------------------------
!
      DO i = 1,nx
        lowi_lut(i) = (i-1)/iskip*iskip + 1
        IF( i == nx) lowi_lut(i) = lowi_lut(i) - iskip
      END DO ! i

      DO j = 1,ny
        lowj_lut(j) = (j-1)/iskip*iskip + 1
        IF( j == ny) lowj_lut(j) = lowj_lut(j) - iskip
      END DO ! i

      DO j=1,ny
        jl = lowj_lut(j)
        jh = jl + iskip
        fracj = FLOAT(j-jl)/FLOAT(iskip)

        DO i=1,nx
          il = lowi_lut(i)
          ih = il + iskip
          fraci = FLOAT(i-il)/FLOAT(iskip)

          z1=t(il,jl,k)
          z2=t(ih,jl,k)
          z3=t(ih,jh,k)
          z4=t(il,jh,k)

          t(i,j,k) =  z1+(z2-z1)*fraci+(z4-z1)*fracj                    &
                    - (z2+z4-z3-z1)*fraci*fracj
        END DO ! i
      END DO ! j

      IF (l_stn_clouds .AND. n_cld_snd >= 1) THEN
!
!-----------------------------------------------------------------------
!
!  Reinterpolate at the station grid points which were skipped
!  previously.
!
!-----------------------------------------------------------------------
!
        DO n1=1,n_cld_snd
          i = i_snd(n1)
          j = j_snd(n1)
          IF(i < 1.OR.i > nx.OR.j < 1.OR.j > ny) GO TO 252
          IF( (i-1) /= (i-1)/iskip*iskip .OR. (j-1) /= (j-1)/iskip*iskip) THEN
            sum=0.
            sumwt=0.
            DO n=nstart,nstop
              ii=iob(n)
              jj=job(n)
              nn=nob(n)
              iii = (i-ii)*(i-ii) + (j-jj)*(j-jj)
              rr = FLOAT (iii)
              rr = spcng*SQRT (rr)
              IF (rr <= range_cld) THEN
!
!-----------------------------------------------------------------------
!
!c                Station obs are being weighted
!
!-----------------------------------------------------------------------
!
                weight = iiilut(i-ii,j-jj) * wt_snd(nn,k)
                sum=weight*cld_snd(nn,k)+sum
                sumwt=sumwt+weight
              END IF
            END DO ! n

            sum = sum + weight_modelfg * cf_modelfg(i,j,k)
            sumwt = sumwt + weight_modelfg

            IF (ABS(sumwt) <= 1.0E-15) THEN
!                  t(i,j,k) = r_missing
!                  istatus = 0
!                  print*,'(00*): ',i,j,k,' tot sum:',sum,
!     :                 '  tot sumwt:',sumwt,'  t:',t(i,j,k)
              t(i,j,k) = 0.01
            ELSE
              t(i,j,k)=sum/sumwt
            END IF

          END IF  ! Stn (i,j) was skipped during the 1st-pass intp.
          252           CONTINUE

        END DO ! n1
      END IF  ! l_stn_clouds = .true. and n_cld_snd >0

    ELSE       ! no 1st guess and obs. is identical to lower lvl
!
!-----------------------------------------------------------------------
!
!c        Obs are identical; Copy analysis from 1 level below
!
!-----------------------------------------------------------------------
!
      IF(cldprt) WRITE(6,51) k,nstart,nstop,nobs
      51        FORMAT(' lvl,nstart,nstop,nobs=',4I5,                   &
                              ' Identical Obs; Copy from 1 lvl down')

      km1 = k - 1

      DO j=1,ny
        DO i=1,nx
          t(i,j,k) = t(i,j,km1)
        END DO ! i
      END DO ! j

    END IF

  END DO ! K

  IF (istatus == 0) THEN
    PRINT*,' WARNING: Distant obs given zero wgt in barnes_r5'
    WRITE(6,*)' Try increasing bias_iii'
  END IF

  istatus = 1
!
  deallocate (fnorm,stat=astat)
  deallocate (iiilut,stat=astat)
  RETURN
END SUBROUTINE barnes_r5





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


SUBROUTINE hgt_to_zk (z,zk,nk,zs_1d,istatus) 2
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  To find the k-position of a given height in the model grid
!  If the height is above the model highest level, the largest
!  level index (nk) will be returned. If the height is below
!  the lowest model level, the lowestlevel index (1) will be
!  returned.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  05/01/96
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  INPUT:
  REAL :: z         ! the height (m) to be coverted to k-index
  INTEGER :: nk     ! number of vertical levels in model grid
  REAL :: zs_1d(nk) ! the physical height of each model level
!
!  OUTPUT:
  INTEGER :: istatus
  REAL :: zk        ! the k-position of the ht. in the model grid
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: k,nkm1
  REAL :: height_above,thickness,frac
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  zk = FLOAT(nk-1)
  height_above = zs_1d(nk-1)

  IF( z > height_above) THEN
    nkm1=nk-1
!        write(6,*)
!        write(6,601) z,zk,nkm1,zs_1d(nkm1)
!601     format(' Warning: above model domain in HGT_TO_ZK'/
!     :        1x,' z=',f9.0,' zk=',f6.2,' nk=',i2,'  zs(nk)=',f9.0)
!        write(6,*) 'Top level index will be assigned to ZK.'
    istatus=0
!
    GO TO 999
  END IF

  DO k = nk-2,1,-1

    IF( height_above >= z .AND. zs_1d(k) <= z) THEN
      thickness = height_above - zs_1d(k)
      frac = (z - zs_1d(k))/thickness
      zk = FLOAT(k) + frac
      istatus=1
      GO TO 999
    END IF

    height_above = zs_1d(k)
  END DO ! K

  zk = 1.0

!      write(6,*)
!      write(6,602) z,zk,zs_1d(1)
!602   format(' Warning: below model level_1 in HGT_TO_ZK'/
!     :      1x,' z=',f9.0,' zk=',f6.2,' k=1  zs(1)=',f9.0)
!      write(6,*) 'Bottom level index will be assigned to ZK.'
  istatus=0

  999   CONTINUE
  RETURN
END SUBROUTINE hgt_to_zk



SUBROUTINE wmix (p,t,dd,w,nl) 1
!
!            07-jun-84       source code from wang
!                            (not on u. of w. source tape)
!
!            dd      - dewpoint depression
  DIMENSION p(*), t(*), dd(*), w(*)
  DO i = 1, nl
    td = t(i) - dd(i)
    IF (td > 253.0) GO TO 110
    es = vpice(td)
    GO TO 120
    110     es = satvap(td)
    120     w(i) = 622.0 * es / p(i)
  END DO
  RETURN
END SUBROUTINE wmix
!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION WSAT                           ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION wsat(p, t),1
!
!         1 4-may-84 added from csu vas code
!                            (not on u. of w. source tape)
!
! jz 4/29/98 to avoid warning message when compile on origin
  DIMENSION p(1),t(1),dd(1),w(1)
  dd(1)  = 0.0
! jz 4/29/98 end
  CALL wmix(p,t,dd,w,1)
  wsat=w(1)
!      wsat=w
!
!         1 4-may-84 return added
!
  RETURN
  END FUNCTION wsat
!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION VPICE                          ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION vpice(temp)
!
!            07-jun-84       source code from wong
!                            (not on u. of w. source tape)
!
  DIMENSION c(6)
!
  REAL*8 c, t
!
  DATA c/ 0.7859063157D+00,  0.3579242320D-01,                          &
       -0.1292820828D-03,  0.5937519208D-06,                            &
        0.4482949133D-09,  0.2176664827D-10/
!
  t = temp - 273.16
  vplog = c(1) + t * (c(2) + t * (c(3) + t * (c(4) + t * (c(5) +        &
               t *  c(6)))))
  vpice = 10.0 ** vplog
  RETURN
  END FUNCTION vpice
!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION SATVAP                         ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION satvap(temp)
!
!            07-jun-84       source code from wang
!                            (not on u. of w. source tape)
!
  DIMENSION c(10)
!
  REAL*8 c, s, t
!
  DATA c/ 0.9999968760D-00, -0.9082695004D-02,                          &
        0.7873616869D-04, -0.6111795727D-06,                            &
        0.4388418740D-08, -0.2988388486D-10,                            &
        0.2187442495D-12, -0.1789232111D-14,                            &
        0.1111201803D-16, -0.3099457145D-19/
  t = temp - 273.16
  s = 6.1078000000D+00 / (c(1) + t * (c(2) + t * (c(3) +                &
                               t * (c(4) + t * (c(5) +                  &
                               t * (c(6) + t * (c(7) +                  &
                               t * (c(8) + t * (c(9) +                  &
                               t * c(10)))))))))) ** 8
  satvap = s
  RETURN
  END FUNCTION satvap


SUBROUTINE precw(p,w,u,np)
!
!            02-may-84       sequence numbers removed
!
  DIMENSION p(*),w(*),u(*)
  DATA f/1961.33/
  w1=w(1)
  p1=p(1)
  s=0.
  u(1)=s
!
!            07-jun-84       correction of tape read error
!  do 100 i=2,n
  DO i = 2, np
    w2=w(i)
    p2=p(i)
    dp=ABS(p2-p1)
    s=s+(w1+w2)*dp/f
    u(i)=s
    w1=w2
    p1=p2
  END DO
  RETURN
END SUBROUTINE precw
!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION VSKINT                         ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION vskint(tbb,jcof,juse,jact)
!
! * obtain sfc-skin-temp from vas brightness temperatures
!     jcof=0 for empirical coefficients
!     jcof=1 for climatological(theoretical) coeffs.
!
!     juse=0 to use 'day'(2-chan) or 'night'(3-chan) eqn.
!        according to solar zenith angle
! *******  currently forced out; will use day only if chan 12 missing.
!     juse=1 to use 'day' eqn.
!     juse=2 to use 'night' eqn.
!
!     jact=0 if nothing done
!     jact=1 if 'day' eqn. actually used
!     jact=2 if 'night' eqn. actually used
  COMMON/nav/vlat,vlon,vzen,szen,il,ie,iras,ipic,itime,jtime,jday
  COMMON/tsurfc/tscc(4,2)
  COMMON/tsurfe/tsce(4,2)
  DIMENSION tsc(4,2),tx(3),ic(3),tbb(*)
  DATA ic/7,8,12/,nx/3/,vmisg/999999./
  DO i=1,nx
    j=ic(i)
    tx(i)=tbb(j)
  END DO
  jact=0
  ts=tx(2)
  IF(tx(1) == vmisg) GO TO 180
  IF(tx(2) == vmisg) GO TO 180
  IF(jcof /= 0) GO TO 110
  DO j=1,2
    DO i=1,4
      tsc(i,j)=tsce(i,j)
    END DO
  END DO
  GO TO 130
  110 DO j=1,2
    DO i=1,4
      tsc(i,j)=tscc(i,j)
    END DO
  END DO
  130 j=juse
  IF(j /= 0) GO TO 150
  j=2                       !changed 1 --> 2
!  if(szen.ge.90.) j=2       !deleted line
  150 IF(j == 1) GO TO 160
  IF(tx(3) == vmisg) j=1
  IF(tx(3) < tx(2)-4.) j=1 !changed tx(2) --> tx(2)-4.
  160 jact=j
  ts=tsc(4,j)
  DO i=1,nx
    ts=ts+tx(i)*tsc(i,j)
  END DO
  IF(ts == 0.) ts=vmisg
  180 vskint=ts
  RETURN
  END FUNCTION vskint
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE PCP_MXR                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE pcp_mxr (nx,ny,nz,t_3d,p_3d ,ref_3d                          & 1
           ,cldpcp_type_3d                                              &
           ,qr_3d,qs_3d,qh_3d,istatus )

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Perform 3D precipitation mixing ratio (in g/kg) analysis using
!  radar reflectivity data. For rain water, using Kessler (1969)
!  formula:
!            qr(g/kg) = a*(rho*arg)**b                  (1)
!
!  Here arg = Z (mm**6/m**3), and dBZ = 10log10 (arg).
!  Coeffcients a=17300.0, and b=7/4.
!  rho represents the air density.
!
!  For snow and hail, using Rogers and Yau (1989) formula:
!
!            qs(g/kg) = c*(rho*arg)**d                  (2)
!
!  where, c=38000.0,  d=2.2
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  06/13/96
!
!  MODIFICATION HISTORY:
!  07/30/97 (J. Zhang)
!           Added precipitation type in the argument list so that
!           mixing ratios of different precip. types can be computed.
!  09/04/97 (J. Zhang)
!            Changed the radar echo thresholds for inserting precip.
!            from radar reflectivities.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  INPUT:
  INTEGER :: nx,ny,nz              ! Model grid size
!
  REAL :: ref_3d(nx,ny,nz)  ! radar reflectivity (dBZ)
  REAL :: t_3d(nx,ny,nz)  ! Temperature (deg. Kelvin)
  REAL :: p_3d(nx,ny,nz)  ! Pressure (Pascal)

  INTEGER*2 cldpcp_type_3d(nx,ny,nz) ! cloud/precip type field
!
!  OUTPUT:
  INTEGER :: istatus
!
  REAL :: qr_3d(nx,ny,nz)     ! rain mixing ratio in (g/kg)
  REAL :: qs_3d(nx,ny,nz)     ! snow/sleet/frz-rain mixing ratio
                            ! in (g/kg)
  REAL :: qh_3d(nx,ny,nz)     ! hail mixing ratio in (g/kg)
!
!  LOCAL:
  REAL :: a,b,c,d                  ! Coef. for Z-qr relation.
  PARAMETER (a=17300.0, b=7.0/4.0)
  PARAMETER (c=38000.0, d=2.2)
  REAL :: rair                     ! Gas constant (J/deg/kg)
  PARAMETER (rair = 287.04)
  REAL :: thresh_ref
  PARAMETER (thresh_ref = 0.0)
  INTEGER :: pcptype
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k, iarg
  REAL :: arg,rhobar,br,dr
  PARAMETER (br=1.0/b, dr=1.0/d)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
  istatus=0
!
!-----------------------------------------------------------------------
!
!  Compute the precip mixing ratio in g/kg from radar reflectivity
!  factor following Kessler (1969) or Rogers and Yau (1989).
!
!-----------------------------------------------------------------------
!
  DO k = 1,nz-1
    DO j = 1,ny-1
      DO i = 1,nx-1
        IF (ref_3d(i,j,k) >= thresh_ref) THEN    ! valid radar refl.
          rhobar = p_3d(i,j,k)/rair/t_3d(i,j,k)
          arg = 10.0**(0.1*ref_3d(i,j,k))
          iarg = cldpcp_type_3d(i,j,k)
          pcptype = iarg/16              ! precip. type

          IF (pcptype == 0) THEN       ! no precip
            PRINT*,'+++ NOTE: radar echo though no precip. +++'
          ELSE IF (pcptype == 1.OR.pcptype == 3) THEN   ! rain or Z R
            qr_3d(i,j,k) = (arg/a)**br/rhobar
          ELSE IF (pcptype == 2) THEN                   ! snow
            qs_3d(i,j,k) = (arg/c)**dr/rhobar
          ELSE IF (pcptype == 4.OR.pcptype == 5) THEN   ! hail or sleet
            qh_3d(i,j,k) = (arg/c)**dr/rhobar
          ELSE                                          ! unknown
            PRINT*,'+++ NOTE: unknown precip type. +++'
          END IF
        ELSE
          qr_3d(i,j,k) = 0.
          qs_3d(i,j,k) = 0.
          qh_3d(i,j,k) = 0.
        END IF
      END DO ! k
    END DO ! i
  END DO ! j
!
!-----------------------------------------------------------------------
!
  istatus = 1
!
  RETURN
END SUBROUTINE pcp_mxr

!wdt-caps sra
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE PCP_MXR_FERRIER             ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE pcp_mxr_ferrier (nx,ny,nz,t_3d,p_3d ,ref_3d                  & 1
           ,cldpcp_type_3d                                              &
           ,qr_3d,qs_3d,qh_3d,istatus )

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Perform 3D precipitation mixing ratio (in g/kg) analysis using
!  radar reflectivity data. For rain water, using Ferrier et al (1995)
!  formulation:
!
!
!     For rain water:
!
!          18
!        10   * 720                              1.75
!  Zer = --------------------------- * (rho * qr)
!          1.75      0.75       1.75
!        pi     * N0r     * rhor
!
!
!     For dry snow (t <= 0 C):
!
!
!          18           2                     0.25
!        10  * 720 * |K|                * rhos
!                       ice                                    1.75
!  Zes = ----------------------------------------- * (rho * qs)     t <= 0 C
!          1.75         2          0.75       2
!        pi        * |K|      * N0s     * rhoi
!                     water
!
!
!     For wet snow (t >= 0 C):
!
!
!              18
!            10   * 720                                 1.75
!  Zes =     ----------------------------   * (rho * qs)            t >  0 C
!              1.75      0.75       1.75
!            pi     * N0s     * rhos
!
!
!     For hail water:
!
!
!          /   18                       \ 0.95
!         /  10   * 720                  \              1.6625
!  Zeh = |   ---------------------------- | * (rho * qh)
!         \    1.75      0.75       1.75 /
!          \ pi     * N0h     * rhoh    /
!
!  Here Zx (mm**6/m**3, x=r,s,h), and dBZ = 10log10 (Zx).
!  rho represents the air density, rhor,rhos,rhoh are the density of
!  rain, snow and hail respectively. Other variables are all constants
!  for this scheme, see below.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Donghai Wang and Eric Kemp)
!  07/20/2000
!
!  MODIFICATION HISTORY:
!
!  11/09/2000 Keith Brewster
!  Moved some parameters with real-valued exponentiation to be
!  computed at runtime due to compiler complaint.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  INPUT:
  INTEGER :: nx,ny,nz              ! Model grid size
!
  REAL :: ref_3d(nx,ny,nz)  ! radar reflectivity (dBZ)
  REAL :: t_3d(nx,ny,nz)  ! Temperature (deg. Kelvin)
  REAL :: p_3d(nx,ny,nz)  ! Pressure (Pascal)

  INTEGER*2 cldpcp_type_3d(nx,ny,nz) ! cloud/precip type field
!
!  OUTPUT:
  INTEGER :: istatus
!
  REAL :: qr_3d(nx,ny,nz)     ! rain mixing ratio in (g/kg)
  REAL :: qs_3d(nx,ny,nz)     ! snow/sleet/frz-rain mixing ratio
                            ! in (g/kg)
  REAL :: qh_3d(nx,ny,nz)     ! hail mixing ratio in (g/kg)
  REAL :: rair                     ! Gas constant (J/deg/kg)
  PARAMETER (rair = 287.04)
  REAL :: thresh_ref
  PARAMETER (thresh_ref = 0.0)
  INTEGER :: pcptype
  REAL :: zerf,zesnegf,zesposf,zehf


  REAL :: a4over7
  PARAMETER (a4over7 = 4.0/7.0)
  REAL :: zehpower
  PARAMETER (zehpower = 1.0/1.6625)
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k, iarg
  REAL :: ze,rho,tc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Intiailize constant factors in the Ze terms for rain, snow and hail,
!  respectively, in Ferrier.
!
!-----------------------------------------------------------------------
!
  istatus=0
  zerf = (1.0E+18 * 720 )                                    &
             / (3.1415926**1.75 * 8.0E+06 ** 0.75 * 1000 ** 1.75 )      &
             / ( 1000. ** 1.75)

  zesnegf = ((1.0E+18 * 720 * 0.176 * 100 ** 0.25 )          &
             /(3.1415926**1.75 * 0.93 * 3.0E+06 ** 0.75 * 917 ** 2))    &
             / ( 1000. ** 1.75)

  zesposf = ((1.0E+18 * 720 )                                &
             / (3.1415926**1.75 * 3.0E+06 ** 0.75 * 100 ** 1.75))       &
             / ( 1000. ** 1.75)

  zehf = (((1.0E+18 * 720 )                                  &
      / (3.1415926**1.75 * 4.0E+04 ** 0.75 * 913 ** 1.75 )) ** 0.95)    &
      / ( 1000. ** 1.6625)

!-----------------------------------------------------------------------
!
!  Compute the precip mixing ratio in g/kg from radar reflectivity
!  factor following Ferrier et al (1995).
!
!-----------------------------------------------------------------------
!
  DO k = 1,nz-1
    DO j = 1,ny-1
      DO i = 1,nx-1
        IF (ref_3d(i,j,k) >= thresh_ref) THEN    ! valid radar refl.
          rho = p_3d(i,j,k)/rair/t_3d(i,j,k)
          ze = 10.0**(0.1*ref_3d(i,j,k))
          iarg = cldpcp_type_3d(i,j,k)
          pcptype = iarg/16              ! precip. type
          tc = t_3d(i,j,k) - 273.15

          IF (pcptype == 0 .OR. rho <= 0.0) THEN       ! no precip
            PRINT*,'+++ NOTE: radar echo though no precip. +++'
          ELSE IF (pcptype == 1.OR.pcptype == 3) THEN   ! rain or Z R
!        print*,'calculating qr...'
            qr_3d(i,j,k) = (ze / zerf ) ** a4over7 / rho
          ELSE IF (pcptype == 2) THEN                   ! snow
            IF (tc <= 0.0) THEN     !dry snow
!          print*,'calculating dry qs...'
              qs_3d(i,j,k) = (ze / zesnegf) ** a4over7 / rho
            ELSE IF (tc > 0.0) THEN     !wet snow
!          print*,'calculating wet qs...'
              qs_3d(i,j,k) = (ze / zesposf) ** a4over7 / rho
              qr_3d(i,j,k) = (ze / zerf ) ** a4over7 / rho
              qr_3d(i,j,k) = qr_3d(i,j,k) - qs_3d(i,j,k)
            END IF
          ELSE IF (pcptype == 4) THEN   ! sleet
!          print*,'calculating wet sleet qs...'
            qs_3d(i,j,k) = (ze / zesposf) ** a4over7 / rho
            qr_3d(i,j,k) = (ze / zerf ) ** a4over7 / rho
            qr_3d(i,j,k) = qr_3d(i,j,k) - qs_3d(i,j,k)
          ELSE IF (pcptype == 5) THEN   ! hail
!      else if (pcptype.eq.4.or.pcptype.eq.5) then   ! hail or sleet
!        print*,'calculating qh...'
            qh_3d(i,j,k) = (ze / zehf) ** zehpower / rho
          ELSE                                          ! unknown
            PRINT*,'+++ NOTE: unknown precip type. +++'
          END IF
        ELSE
          qr_3d(i,j,k) = 0.
          qs_3d(i,j,k) = 0.
          qh_3d(i,j,k) = 0.
        END IF
      END DO ! k
    END DO ! i
  END DO ! j
  PRINT*,'Finish Ferrier ...'
!
!-----------------------------------------------------------------------
!
  istatus = 1
!
  RETURN
END SUBROUTINE pcp_mxr_ferrier
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE REFMOSAIC                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!

SUBROUTINE refmosaic(nradfil,nx,ny,nz,mx_rad,                           & 1,1
           xs,ys,zs,radfname,ref_mos_3d,rhinf,rvinf,                    &
           refl,tem1,tem2,istatus)
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  This routine patches together reflectivity fields observed from
!  multiple radars to create one 3D radar reflectivity field.
!  It also fills the data gaps between the radar beams and elevations).
!  These gaps are possible because that the real-time radar data
!  processing only put data to the grid point nearest to the radar gate.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  06/04/96.
!
!  MODIFICATION HISTORY:
!  11/01/01 (Keith Brewster)  
!  Created refmosaic from rad_patch_fill saving lots of memory when
!  there are more than a couple of radars.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
  INTEGER :: nradfil              ! # of radars currently available
  INTEGER :: nx,ny,nz             ! model grid size
  INTEGER :: mx_rad
  REAL :: xs(nx)                  ! x-coor at scalar points
  REAL :: ys(ny)                  ! y-coor at scalar points
  REAL :: zs(nx,ny,nz)            ! z-coor at scalar points
  CHARACTER (LEN=132) radfname(mx_rad)
  REAL :: rhinf
  REAL :: rvinf
!
!  OUTPUT:
!
  REAL :: ref_mos_3d (nx,ny,nz) ! combined radar refl. field
  INTEGER :: istatus
!
!  Temporary working array
!
  REAL :: refl(nx,ny,nz)
  REAL :: tem1(nx,ny,nz)
  REAL :: tem2(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,irad
  INTEGER :: im1,ip1,jm1,jp1,km1,kp1
  INTEGER :: isrcrad,istat_rad
  REAL    :: latrad,lonrad,elvrad
  REAL    :: x_m1,x_p1,y_m1,y_p1,hgt_km1,hgt_kp1,frac
  REAL    :: refl_ip1,refl_im1,refl_jm1,refl_jp1,refl_km1,refl_kp1
  LOGICAL :: got_im1,got_ip1,got_jm1,got_jp1,got_km1,got_kp1
  INTEGER :: strhopt_rad,mapproj_rad
  REAL    :: dx_rad,dy_rad,dz_rad,dzmin_rad
  REAL    :: ctrlat_rad,ctrlon_rad,tlat1_rad,tlat2_rad,tlon_rad
  REAL    :: sclfct_rad
  CHARACTER (LEN=5) stnrad
  CHARACTER (LEN=128) warn_string 
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
   INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  istatus=0

  warn_string = 'WARNING: Radar data are not used ' //                  &
                   'due to grid or map projection inconsistencies'

  WRITE(6,'(a)')                                                        &
    ' Mosaiking reflectivity from multiple radar observations'
!
!-----------------------------------------------------------------------
!
!  Combining the multiple radar obs. to produce one set of
!  3D reflectivity field.
!
!-----------------------------------------------------------------------
!
  DO irad = 1,nradfil
    CALL read1rad(nx,ny,nz,radfname(irad),                              &
                  strhopt_rad,mapproj_rad,dx_rad,dy_rad,dz_rad,         &
                  dzmin_rad,ctrlat_rad,ctrlon_rad,                      &
                  tlat1_rad,tlat2_rad,tlon_rad,sclfct_rad,              &
                  isrcrad,stnrad,latrad,lonrad,elvrad,                  &
                  tem1,refl,tem2,tem2,                                  &
                  istat_rad)
!
!-----------------------------------------------------------------------
!
!  For successful read of radar data.
!
!  Check IF the radar data are remapped on the same grid as
!  the one we are using for ADAS.
!
!-----------------------------------------------------------------------
!
    IF (istat_rad == 1 ) THEN
      WRITE(6,'(a,a)') ' Successfully read radar datafile ',radfname(irad)
      IF (ABS(strhopt_rad-strhopt) > 1.0E-4 .OR.                        &
              mapproj_rad /= mapproj) THEN
        WRITE(6,*) warn_string
        WRITE(6,*)                                                      &
          'RADAR: strhopt=',strhopt_rad, '  mapproj=',mapproj_rad
        WRITE(6,*) 'ADAS: strhopt=',strhopt, '  mapproj=',mapproj
        CYCLE
      END IF
!
      IF (ABS(dx_rad-dx) > 1.0E-4 .OR. ABS(dy_rad-dy) > 1.0E-4 .OR.     &
          ABS(dz_rad-dz) > 1.0E-4 .OR.                                  &
          ABS(dzmin_rad-dzmin) > 1.0E-4) THEN
        WRITE(6,*) warn_string
        WRITE(6,*) 'RADAR: dx=',dx_rad, ' dy=',dy_rad,' dz=',dz_rad
        WRITE(6,*) 'ADAS: dx=',dx, ' dy=',dy,' dz=',dz
        CYCLE
      END IF
!
      IF (mapproj /= 0) THEN
        IF (ABS(ctrlat_rad-ctrlat) > 1.0E-4 .OR.                        &
              ABS(ctrlon_rad-ctrlon) > 1.0E-4) THEN
          WRITE(6,*) warn_string
          WRITE(6,*) 'RADAR: ctrlat=',ctrlat_rad,' ctrlon=',ctrlon_rad
          WRITE(6,*) 'ADAS: ctrlat=',ctrlat, ' ctrlon=',ctrlon
          CYCLE
        END IF
!
        IF(ABS(sclfct_rad-sclfct) > 1.0E-4) THEN
          WRITE(6,*) warn_string
          WRITE(6,*) 'RADAR: sclfct=',sclfct_rad
          WRITE(6,*) 'ADAS: sclfct=',sclfct
          CYCLE
        END IF
      END IF
!
      IF (mapproj == 1 .OR. mapproj == 3) THEN
        IF (ABS(tlat1_rad-trulat1) > 1.0E-4 .OR.                        &
              ABS(tlon_rad-trulon) > 1.0E-4) THEN
          WRITE(6,*) warn_string
          WRITE(6,*) 'RADAR: trulat1=',tlat1_rad, ' trulon=',tlon_rad
          WRITE(6,*) 'ADAS: trulat1=',trulat1, ' trulon=',trulon
          CYCLE
        END IF
      ELSE IF (mapproj == 2) THEN
        IF (ABS(tlat1_rad-trulat1) > 1.0E-4 .OR.                        &
            ABS(tlat2_rad-trulat2) > 1.0E-4 &
            .OR. ABS(tlon_rad-trulon) > 1.0E-4) THEN
          WRITE(6,*) warn_string
          WRITE(6,*) 'RADAR: trulat1=',tlat1_rad,                       &
                   ' trulat2=',tlat2_rad,' trulon=',tlon_rad
          WRITE(6,*) 'ADAS: trulat1=',trulat1,                          &
                   ' trulat2=',trulat2,' trulon=',trulon
          CYCLE
        END IF
      END IF
!
!-----------------------------------------------------------------------
!
!     For now, mosaic by taking the maximum value
!
!-----------------------------------------------------------------------
!
      WRITE(6,'(a,a,a)') ' Bringing radar ',stnrad,' into mosaic '

! OpenMP changed loop order to j,k,i:
!$OMP PARALLEL DO PRIVATE (j,k,i,irad)
      DO j=1,ny
        DO k=1,nz-1
          DO i=1,nx
            ref_mos_3d(i,j,k) = MAX(ref_mos_3d(i,j,k),refl(i,j,k))
          END DO
        END DO
      END DO
    END IF ! read ok
  END DO  !irad

!
!-----------------------------------------------------------------------
!
!  Vertical interpolation to fill the gaps between different
!  elevations.
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(a)') ' Vertically interpolating radar refl.'
! OpenMP changed loop order to j,i,k:
!$OMP PARALLEL DO PRIVATE(j,i,k,kp1,got_kp1,refl_kp1,hgt_kp1,km1,got_km1,refl_km1,hgt_km1,frac)
  DO j=1,ny
    DO i=1,nx
      DO k=2,nz-2
        tem1(i,j,k) = ref_mos_3d(i,j,k)
        IF(ref_mos_3d(i,j,k) < -10.) THEN
          got_kp1=.FALSE.
          DO kp1=k+1,nz-1
            IF (zs(i,j,kp1) > (zs(i,j,k)+rvinf)) EXIT
            IF (ref_mos_3d(i,j,kp1) >= -10.) THEN
              refl_kp1 = ref_mos_3d(i,j,kp1)
              hgt_kp1 = zs(i,j,kp1)
              got_kp1=.TRUE.
              EXIT
            END IF
          END DO

          got_km1=.FALSE.
          DO km1=k-1,2,-1
            IF(zs(i,j,km1) < (zs(i,j,k)-rvinf)) EXIT
            IF(ref_mos_3d(i,j,km1) >= -10) THEN
              refl_km1 = ref_mos_3d(i,j,km1)
              hgt_km1 = zs(i,j,km1)
              got_km1=.TRUE.
              EXIT
            END IF
          END DO

          IF(got_km1 .AND. got_kp1 .AND.                             &
             (hgt_kp1-hgt_km1) < rvinf ) THEN
            frac = (zs(i,j,k)-hgt_km1)/(hgt_kp1-hgt_km1)
            tem1(i,j,k) = refl_kp1*frac + refl_km1*(1.-frac)
          END IF
        END IF
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Horizontal interpolation to further fill gaps between
!  the beams.
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(a)') ' Horizontally interpolating radar refl.'

  DO k=1,nz-1
    DO i=2,nx-2
      DO j=2,ny-2
        ref_mos_3d(i,j,k)=tem1(i,j,k)
        IF(tem1(i,j,k) < -10.) THEN 
          got_ip1=.FALSE.
          DO ip1=i+1,nx-1
            IF (xs(ip1) > (xs(i)+rhinf)) EXIT
            IF(tem1(ip1,j,k) >= -10.) THEN
              refl_ip1=tem1(ip1,j,k)
              x_p1=xs(ip1)
              got_ip1=.TRUE.
              EXIT
            END IF
          END DO

          got_im1=.FALSE.
          DO im1=i-1,1,-1
            IF (xs(im1) < (xs(i)-rhinf)) EXIT
            IF(tem1(im1,j,k) >= -10.) THEN
              refl_im1=tem1(im1,j,k)
              x_m1=xs(im1)
              got_im1=.TRUE.
              EXIT
            END IF
          END DO

          got_jp1=.FALSE.
          DO jp1=j+1,ny-1
            IF (ys(jp1) > (ys(j)+rhinf)) EXIT
            IF(tem1(i,jp1,k) >= -10.) THEN
              refl_jp1=tem1(i,jp1,k)
              y_p1=ys(jp1)
              got_jp1=.TRUE.
              EXIT
            END IF
          END DO

          got_jm1=.FALSE.
          DO jm1=j-1,1,-1
            IF (ys(jm1) < (ys(j)-rhinf)) EXIT
            IF(tem1(i,jm1,k) >= -10.) THEN
              refl_jm1=tem1(i,jm1,k)
              y_m1=ys(jm1)
              got_jm1=.TRUE.
              EXIT
            END IF
          END DO

          IF(got_im1 .AND. got_ip1 .AND. (x_p1-x_m1) < rhinf .AND.  &
             got_jm1 .AND. got_jp1 .AND. (y_p1-y_m1) < rhinf ) THEN
            frac = (xs(i)-x_m1)/(x_p1-x_m1)
            ref_mos_3d(i,j,k) = refl_ip1*frac + refl_im1*(1-frac)
            frac = (ys(j)-y_m1)/(y_p1-y_m1)
            ref_mos_3d(i,j,k) = (refl_jp1*frac + refl_jm1*(1-frac)    &
                                     + ref_mos_3d(i,j,k))*0.5
          ELSE IF (got_im1 .AND. got_ip1 .AND.                      &
                   (x_p1-x_m1) < rhinf) THEN
            frac = (xs(i)-x_m1)/(x_p1-x_m1)
            ref_mos_3d(i,j,k) = refl_ip1*frac + refl_im1*(1-frac)
          ELSE IF (got_jm1 .AND. got_jp1 .AND.                      &
                   (y_p1-y_m1) < rhinf) THEN
            frac = (ys(j)-y_m1)/(y_p1-y_m1)
            ref_mos_3d(i,j,k) = refl_jp1*frac + refl_jm1*(1-frac)
          END IF
        END IF
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
! Done.
!
!-----------------------------------------------------------------------
!
  istatus=1

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


SUBROUTINE readrad(nx,ny,nz,isrcrad,stnrad                              &,4
           ,latrad,lonrad,elvrad                                        &
           ,gridvel,gridref,gridnyq,gridtim                             &
           ,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Reads radar data remapped on the ARPS grid.
!  This routine requires the remapping to occur on the same grid
!  as the analysis.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Jian Zhang
!  05/1996  Read the remapped radar data which was written by the
!           corresponding output routine "wrtrad" in remaplib.f.
!
!  MODIFICATION HISTORY:
!  03/19/97  J. Zhang
!            Added a line of error message when there is trouble
!            reading a radar file.
!  04/03/97  J. Zhang
!            Added the option of reading the data file created
!            from "WTRADCOL".  Added output for the remapping
!            parameters in the radar file (e.g., strhopt,mapproj,
!            dx,dy,dz,dzmin,ctrlat,ctrlon,tlat1,tlat2,tlon,scale)
!  04/07/97  J. Zhang
!            Added  the QC for the case when i,j,k outside the model
!            domain
!  04/09/97  J. Zhang
!            Added the Initializations for gridref, girdvel...
!  04/11/97  J. Zhang
!            Include dims.inc for nx,ny,nz
!  04/14/97  J. Zhang
!            Added message output for the case when actual # of
!            radar files exceeds the maximum allowed number in the
!            ADAS include file.  When that happens, the program will
!            stop.
!  03/31/98  J. Zhang
!            Deleted the option for reading the radar data file
!            created from "WRTRAD".
!
!-----------------------------------------------------------------------
!
!  INCLUDE:  (from dims.inc and adas.inc)
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!    mx_rad     maximum number of radars
!
!    nradfil    number of radar files
!    radfname   file name for radar datasets
!
!  OUTPUT:
!
!    isrcrad  index of radar source
!    stnrad   radar site name    character*4
!    latrad   latitude of radar  (degrees N)
!    lonrad   longitude of radar (degrees E)
!    elvrad   elevation of feed horn of radar (m MSL)
!
!    gridvel  radial velocity on ARPS grid
!    gridref  reflectivity on ARPS grid
!    gridnyq  nyquist velocity on ARPS grid
!    gridtim  observation time at ARPS grid
!
!    istatus  status indicator
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
  INCLUDE 'adas.inc'     ! ADAS parameters
!
!-----------------------------------------------------------------------
!
!  INPUT:
  INTEGER :: nx,ny,nz    ! the ARPS grid size
!
!  INCLUDE:   (from adas.inc)
!  integer mx_rad
!  INPUT:   (from namelist in the .input file)
!  integer nradfil
!  character*80 radfname(mx_rad)
!
!  LOCAL:
  REAL :: readk(nz)
  REAL :: readhgt(nz)
  REAL :: readref(nz)
  REAL :: readvel(nz)
  REAL :: readnyq(nz)
  REAL :: readtim(nz)
!
  INTEGER :: kntref(nz)
  INTEGER :: kntvel(nz)
  INTEGER :: iradvr
  INTEGER :: nradvr
!
!  OUTPUT:
  INTEGER :: istatus
!
!  OUTPUT:  ARPS radar arrays
  REAL :: gridvel(nx,ny,nz,mx_rad)
  REAL :: gridref(nx,ny,nz,mx_rad)
  REAL :: gridnyq(nx,ny,nz,mx_rad)
  REAL :: gridtim(nx,ny,nz,mx_rad)
!
!  OUTPUT:  Radar site variables
  INTEGER :: isrcrad(0:mx_rad)
  CHARACTER (LEN=5) :: stnrad(mx_rad)
  REAL :: latrad(mx_rad)
  REAL :: lonrad(mx_rad)
  REAL :: elvrad(mx_rad)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=4) :: stn
  CHARACTER (LEN=6) :: runname
  CHARACTER (LEN=132) :: fname
  INTEGER :: ireftim,itime,vcpnum,idummy
  INTEGER :: hdmpfmt,strhopt,mapprin
  INTEGER :: nchanl,ierr
  INTEGER :: iyr, imon, idy, ihr, imin, isec
  INTEGER :: i,j,k,krad,kk,ipt,klev

  REAL :: dxin,dyin,dzin,dzminin,ctrlatin
  REAL :: ctrlonin,tlat1in,tlat2in,tlonin,scalin,rdummy
  REAL :: xrd,yrd,gridlat,gridlon,elev
!
!-----------------------------------------------------------------------
!
!  Common block that stores remapping parameters for the radar
!  data file.
!
!-----------------------------------------------------------------------
!
  COMMON/remapfactrs_rad/strhopt,mapprin
  COMMON/remapfactrs_rad2/dxin,dyin,dzin,dzminin,                       &
           ctrlatin,ctrlonin,tlat1in,tlat2in,tlonin,scalin
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
  istatus=0
  IF (nradfil > mx_rad) THEN
    WRITE(6,'(a,i3,a,i3/a)')                                            &
        ' ERROR: nradfil ',nradfil,' exceeds mx_rad dimension',         &
        mx_rad,' please increase MX_RAD in the .inc file'
    PRINT*,' ABORTING from READRAD......'
    STOP
  END IF
  IF(nradfil < 1) THEN
    WRITE(6,*) 'No radar data available. Returning from READRAD...'
    RETURN
  END IF
!
!-----------------------------------------------------------------------
!
!  Initializations
!
!-----------------------------------------------------------------------
!
! OpenMP changed loop order to j,krad,k,i:
!$OMP PARALLEL DO PRIVATE(j,krad,k,i)
  DO j=1,ny
    DO krad = 1, mx_rad
      DO k=1,nz
        DO i=1,nx
          gridref(i,j,k,krad)=-9999.
          gridvel(i,j,k,krad)=-9999.
          gridnyq(i,j,k,krad)=-9999.
          gridtim(i,j,k,krad)=-9999.
        END DO
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Loop through all radars
!
!-----------------------------------------------------------------------
!
  DO krad = 1, nradfil

    fname=radfname(krad)
    CALL asnctl ('NEWLOCAL', 1, ierr)
    CALL asnfile(fname, '-F f77 -N ieee', ierr)

    CALL getunit( nchanl )
    OPEN(UNIT=nchanl,FILE=trim(fname),ERR=399,                          &
         FORM='unformatted',STATUS='old')
!
!-----------------------------------------------------------------------
!
!  Read radar description variables
!
!-----------------------------------------------------------------------
!
    istatus=1
    isrcrad(krad)=1
!
    READ(nchanl) stn
    stnrad(krad)=stn
    READ(nchanl) ireftim,itime,vcpnum,idummy,idummy,                    &
               idummy,idummy,idummy,idummy,idummy
!
    CALL abss2ctim(itime, iyr, imon, idy, ihr, imin, isec )
    iyr=MOD(iyr,100)
    WRITE(6,'(/a,i2.2,a,i2.2,a,i2.2,1X,i2.2,a,i2.2,a)')                 &
       'Reading remapped raw radar data for: ',                         &
       imon,'/',idy,'/',iyr,ihr,':',imin,' UTC'
!
    READ(nchanl) runname
    READ(nchanl) hdmpfmt,strhopt,mapprin,idummy,idummy,                 &
               idummy,idummy,idummy,idummy,idummy

    READ(nchanl) dxin,dyin,dzin,dzminin,ctrlatin,                       &
           ctrlonin,tlat1in,tlat2in,tlonin,scalin,                      &
           latrad(krad),lonrad(krad),elvrad(krad),                      &
           rdummy,rdummy
!
    DO k=1,nz
      kntref(k) = 0
      kntvel(k) = 0
    END DO

    READ(nchanl) nradvr,iradvr
!
    DO ipt=1,(nx*ny)

      READ(nchanl,END=51) i,j,xrd,yrd,                                  &
                     gridlat,gridlon,elev,klev
      READ(nchanl,END=52) (readk(kk),kk=1,klev)
      READ(nchanl,END=52) (readhgt(kk),kk=1,klev)
      READ(nchanl,END=52) (readref(kk),kk=1,klev)
      READ(nchanl,END=52) (readvel(kk),kk=1,klev)
      READ(nchanl,END=52) (readnyq(kk),kk=1,klev)
      READ(nchanl,END=52) (readtim(kk),kk=1,klev)

      IF(i <= nx.AND.i >= 1 .AND. j <= ny.AND.j >= 1) THEN
        DO kk=1,klev
          k=nint(readk(kk))
          IF(k <= nz.AND.k >= 1) THEN
            gridref(i,j,k,krad)=readref(kk)
            gridvel(i,j,k,krad)=readvel(kk)
            gridnyq(i,j,k,krad)=readnyq(kk)
            gridtim(i,j,k,krad)=readtim(kk)
            IF (gridref(i,j,k,krad) > -200. .AND. gridref(i,j,k,krad) < 200.) &
                kntref(k)=kntref(k)+1
            IF (gridvel(i,j,k,krad) > -200. .AND. gridvel(i,j,k,krad) < 200.) &
                kntvel(k)=kntvel(k)+1
          END IF  ! 1 < k < nz
        END DO  ! kk = 1, klev
      END IF  ! 1 < i < nx  & 1 < j < ny

    END DO  ! ipt = 1, nx*ny

    51      CONTINUE
    ipt=ipt-1
    WRITE(6,'(a,i6,a)') ' End of file reached after reading',           &
                       ipt,' columns'
    GO TO 55
    52      CONTINUE
    WRITE(6,'(a,i6,a)') ' End of file reached while reading',           &
                       ipt,' column'
    55      CONTINUE
!
!-----------------------------------------------------------------------
!
!  Write statistics
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(a)') '  k   n ref    n vel'
    DO k=1,nz
      WRITE(6,'(i3,2i10)') k,kntref(k),kntvel(k)
    END DO
!
    CLOSE(nchanl)
    CALL retunit( nchanl )

    GO TO 400

    399     CONTINUE
    PRINT*,'Error reading the radar file:',fname

    400     CONTINUE
  END DO  ! KRAD = 1, nradfil

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

SUBROUTINE read1rad(nx,ny,nz,radfname,                                  & 1,35
                  strhopt_rad,mapproj_rad,dx_rad,dy_rad,dz_rad,         &
                  dzmin_rad,ctrlat_rad,ctrlon_rad,                      &
                  tlat1_rad,tlat2_rad,tlon_rad,sclfct_rad,              &
                  isrcrad,stnrad,latrad,lonrad,elvrad,                  &
                  gridvel,gridref,gridnyq,gridtim,                      &
                  istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Reads radar data remapped on the ARPS grid.
!  This routine requires the remapping to occur on the same grid
!  as the analysis.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Jian Zhang
!  05/1996  Read the remapped radar data which was written by the
!           corresponding output routine "wrtrad" in remaplib.f.
!
!  MODIFICATION HISTORY:
!  03/19/97  J. Zhang
!            Added a line of error message when there is trouble
!            reading a radar file.
!  04/03/97  J. Zhang
!            Added the option of reading the data file created
!            from "WTRADCOL".  Added output for the remapping
!            parameters in the radar file (e.g., strhopt,mapproj,
!            dx,dy,dz,dzmin,ctrlat,ctrlon,tlat1,tlat2,tlon,scale)
!  04/07/97  J. Zhang
!            Added  the QC for the case when i,j,k outside the model
!            domain
!  04/09/97  J. Zhang
!            Added the Initializations for gridref, girdvel...
!  04/11/97  J. Zhang
!            Include dims.inc for nx,ny,nz
!  04/14/97  J. Zhang
!            Added message output for the case when actual # of
!            radar files exceeds the maximum allowed number in the
!            ADAS include file.  When that happens, the program will
!            stop.
!  03/31/98  J. Zhang
!            Deleted the option for reading the radar data file
!            created from "WRTRAD".
!  11/01/01  K. Brewster
!            Modified readrad to read1rad to read only one radar file 
!            at a time to save memory.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!    radfname   file name for radar datasets
!
!  OUTPUT:
!
!    isrcrad  index of radar source
!    stnrad   radar site name    character*4
!    latrad   latitude of radar  (degrees N)
!    lonrad   longitude of radar (degrees E)
!    elvrad   elevation of feed horn of radar (m MSL)
!
!    gridvel  radial velocity on ARPS grid
!    gridref  reflectivity on ARPS grid
!    gridnyq  nyquist velocity on ARPS grid
!    gridtim  observation time at ARPS grid
!
!    istatus  status indicator
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz    ! the ARPS grid size
  CHARACTER*132 radfname
  INTEGER :: strhopt_rad,mapproj_rad
  REAL    :: dx_rad,dy_rad,dz_rad,dzmin_rad
  REAL    :: ctrlat_rad,ctrlon_rad,tlat1_rad,tlat2_rad,tlon_rad
  REAL    :: sclfct_rad
!
!  OUTPUT:  Radar site variables
!
  INTEGER :: isrcrad
  CHARACTER (LEN=5) :: stnrad
  REAL :: latrad
  REAL :: lonrad
  REAL :: elvrad
!
!  OUTPUT:  ARPS radar arrays
!
  REAL :: gridvel(nx,ny,nz)
  REAL :: gridref(nx,ny,nz)
  REAL :: gridnyq(nx,ny,nz)
  REAL :: gridtim(nx,ny,nz)
  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: readk(nz)
  REAL :: readhgt(nz)
  REAL :: readref(nz)
  REAL :: readvel(nz)
  REAL :: readnyq(nz)
  REAL :: readtim(nz)
!
  INTEGER :: kntref(nz)
  INTEGER :: kntvel(nz)
!
  INTEGER :: mxradvr,nradvr
  PARAMETER(mxradvr=10)
  INTEGER :: iradvr(mxradvr)

  CHARACTER (LEN=4) :: stn
  CHARACTER (LEN=6) :: runname
  INTEGER :: ireftim,itime,vcpnum,idummy
  INTEGER :: hdmpfmt,strhopt,mapprin
  INTEGER :: nchanl,ierr
  INTEGER :: iyr, imon, idy, ihr, imin, isec
  INTEGER :: i,j,k,krad,kk,ipt,klev

  REAL :: xrd,yrd,gridlat,gridlon,elev,rdummy
  INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmp(:,:,:) ! Temporary array
  REAL, ALLOCATABLE :: hmax(:), hmin(:) ! Temporary array
  INTEGER :: lens,dmpfmt,sd_id,isource
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
  istatus=0
!
!-----------------------------------------------------------------------
!
!  Initializations
!
!-----------------------------------------------------------------------
!
! OpenMP changed loop order to j,k,i:
!$OMP PARALLEL DO PRIVATE(j,k,i)
  DO j=1,ny
    DO k=1,nz
      DO i=1,nx
        gridref(i,j,k)=-9999.
        gridvel(i,j,k)=-9999.
        gridnyq(i,j,k)=-9999.
        gridtim(i,j,k)=-9999.
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
! Open file
!
!-----------------------------------------------------------------------
!
  CALL asnctl ('NEWLOCAL', 1, ierr)
  CALL asnfile(radfname, '-F f77 -N ieee', ierr)

  lens=LEN(trim(radfname))
  IF(radfname(lens-3:lens)=='hdf4')THEN
    dmpfmt=2
  ELSE
    dmpfmt=1
  ENDIF
  print*,'radfname==',radfname,'dmpfmt==',dmpfmt

  IF(dmpfmt==1)THEN
  CALL getunit( nchanl )
  OPEN(UNIT=nchanl,FILE=trim(radfname),ERR=399,                         &
         FORM='unformatted',STATUS='old')
!
!-----------------------------------------------------------------------
!
!  Read radar description variables
!
!-----------------------------------------------------------------------
!
  istatus=1
  isrcrad=1
!
  READ(nchanl) stn
  stnrad=stn
  READ(nchanl) ireftim,itime,vcpnum,idummy,idummy,                      &
               idummy,idummy,idummy,idummy,idummy
!
  CALL abss2ctim(itime, iyr, imon, idy, ihr, imin, isec )
  iyr=MOD(iyr,100)
  WRITE(6,'(/a,i2.2,a,i2.2,a,i2.2,1X,i2.2,a,i2.2,a)')                   &
       'Reading remapped raw radar data for: ',                         &
       imon,'/',idy,'/',iyr,ihr,':',imin,' UTC'
!
  READ(nchanl) runname
  READ(nchanl) hdmpfmt,strhopt_rad,mapproj_rad,idummy,idummy,           &
               idummy,idummy,idummy,idummy,idummy

  READ(nchanl) dx_rad,dy_rad,dz_rad,dzmin_rad,ctrlat_rad,               &
           ctrlon_rad,tlat1_rad,tlat2_rad,tlon_rad,sclfct_rad,          &
           latrad,lonrad,elvrad,                                        &
           rdummy,rdummy
!
  DO k=1,nz
    kntref(k) = 0
    kntvel(k) = 0
  END DO

  READ(nchanl) nradvr,iradvr
!
  DO ipt=1,(nx*ny)

    READ(nchanl,END=51) i,j,xrd,yrd,                                    &
                   gridlat,gridlon,elev,klev
    READ(nchanl,END=52) (readk(kk),kk=1,klev)
    READ(nchanl,END=52) (readhgt(kk),kk=1,klev)
    READ(nchanl,END=52) (readref(kk),kk=1,klev)
    READ(nchanl,END=52) (readvel(kk),kk=1,klev)
    READ(nchanl,END=52) (readnyq(kk),kk=1,klev)
    READ(nchanl,END=52) (readtim(kk),kk=1,klev)

    IF(i <= nx.AND.i >= 1 .AND. j <= ny.AND.j >= 1) THEN
      DO kk=1,klev
        k=nint(readk(kk))
        IF(k <= nz.AND.k >= 1) THEN
          gridref(i,j,k)=readref(kk)
          gridvel(i,j,k)=readvel(kk)
          gridnyq(i,j,k)=readnyq(kk)
          gridtim(i,j,k)=readtim(kk)
          IF (gridref(i,j,k) > -200. .AND. gridref(i,j,k) < 200.)      &
              kntref(k)=kntref(k)+1
          IF (gridvel(i,j,k) > -200. .AND. gridvel(i,j,k) < 200.)      &
              kntvel(k)=kntvel(k)+1
        END IF  ! 1 < k < nz
      END DO  ! kk = 1, klev
    END IF  ! 1 < i < nx  & 1 < j < ny

  END DO  ! ipt = 1, nx*ny

  51      CONTINUE
  ipt=ipt-1
    WRITE(6,'(a,i6,a)') ' End of file reached after reading',          &
                       ipt,' columns'
    GO TO 55
    52      CONTINUE
    WRITE(6,'(a,i6,a)') ' End of file reached while reading',          &
                       ipt,' column'
    55      CONTINUE
!
!-----------------------------------------------------------------------
!
!  Write statistics
!
!-----------------------------------------------------------------------
!
    print*,'finishing reading '
    WRITE(6,'(a)') '  k   n ref    n vel'
    DO k=1,nz
      WRITE(6,'(i3,2i10)') k,kntref(k),kntvel(k)
    END DO
!
    CLOSE(nchanl)
    CALL retunit( nchanl )

    GO TO 400

    399     CONTINUE
    PRINT*,'Error reading the radar file:',radfname

    400     CONTINUE
  ELSE !HDF4 file 

    ALLOCATE (itmp(nx,ny,nz),stat=istatus)
    IF (istatus /= 0) THEN
      WRITE (6,*) "HDFREAD: ERROR allocating itmp, returning"
      STOP
    END IF
    ALLOCATE (hmax(nz),stat=istatus)
    IF (istatus /= 0) THEN
      WRITE (6,*) "HDFREAD: ERROR allocating hmax, returning"
      STOP
    END IF
    ALLOCATE (hmin(nz),stat=istatus)
    IF (istatus /= 0) THEN
      WRITE (6,*) "HDFREAD: ERROR allocating hmin, returning"
      STOP
    END IF

    CALL hdfopen(trim(radfname), 1, sd_id)
    IF (sd_id < 0) THEN
    WRITE (6,*) "RDRADCOL: ERROR opening ",                           &
                trim(radfname)," for reading."
    istatus = 1
    STOP
    END IF
    CALL hdfrdc(sd_id, 4, 'radid', stn, istatus)
    stnrad=stn
    CALL hdfrdi(sd_id, 'ireftim', ireftim, istatus)
    CALL hdfrdi(sd_id, 'itime', itime, istatus)
    CALL hdfrdi(sd_id, 'vcpnum', vcpnum, istatus)
    CALL hdfrdi(sd_id, 'isource', isource, istatus)
    CALL hdfrdc(sd_id, 4, 'runname', runname, istatus)
    CALL hdfrdi(sd_id, 'hdmpfmt', hdmpfmt, istatus)
    CALL hdfrdi(sd_id, 'strhopt', strhopt_rad, istatus)
    CALL hdfrdi(sd_id, 'mapproj', mapproj_rad, istatus)
    CALL hdfrdr(sd_id, 'dx', dx_rad, istatus)
    CALL hdfrdr(sd_id, 'dy', dy_rad, istatus)
    CALL hdfrdr(sd_id, 'dz', dz_rad, istatus)
    CALL hdfrdr(sd_id, 'dzmin', dzmin_rad, istatus)
    CALL hdfrdr(sd_id, 'ctrlat', ctrlat_rad, istatus)
    CALL hdfrdr(sd_id, 'ctrlon', ctrlon_rad, istatus)
    CALL hdfrdr(sd_id, 'trulat1', tlat1_rad, istatus)
    CALL hdfrdr(sd_id, 'trulat2', tlat2_rad, istatus)
    CALL hdfrdr(sd_id, 'trulon', tlon_rad, istatus)
    CALL hdfrdr(sd_id, 'sclfct', sclfct_rad, istatus)
    CALL hdfrdr(sd_id, 'latrad', latrad, istatus)
    CALL hdfrdr(sd_id, 'lonrad', lonrad, istatus)
    CALL hdfrdr(sd_id, 'elvrad', elvrad, istatus)
    CALL hdfrdi(sd_id, 'nradvr', nradvr, istatus)
    CALL hdfrd1d(sd_id,'iradvr', mxradvr,iradvr,'')
    PRINT *, ' Got nradvr,iradvr: ',nradvr,iradvr
    CALL abss2ctim(itime, iyr, imon, idy, ihr, imin, isec )
    iyr=MOD(iyr,100)
    WRITE(6,'(/a,i2.2,a,i2.2,a,i2.2,1X,i2.2,a,i2.2,a)')                 &
       'Reading remapped raw radar data for: ',                         &
       imon,'/',idy,'/',iyr,ihr,':',imin,' UTC'

    CALL hdfrd3d(sd_id,"gridref",nx,ny,nz,gridref,istatus,itmp,hmax,hmin)
    IF (istatus /= 0) GO TO 115
    CALL hdfrd3d(sd_id,"gridvel",nx,ny,nz,gridvel,istatus,itmp,hmax,hmin)
    IF (istatus /= 0) GO TO 115
    CALL hdfrd3d(sd_id,"gridnyq",nx,ny,nz,gridnyq,istatus,itmp,hmax,hmin)
    IF (istatus /= 0) GO TO 115
    CALL hdfrd3d(sd_id,"gridtim",nx,ny,nz,gridtim,istatus,itmp,hmax,hmin)
    IF (istatus /= 0) GO TO 115
    CALL hdfclose(sd_id,istatus)
    DO i=1,nx
    DO j=1,ny
    DO k=1,nz
      IF (gridref(i,j,k) > -200. .AND. gridref(i,j,k) < 200.)     &
        kntref(k)=kntref(k)+1
      IF (gridvel(i,j,k) > -200. .AND. gridvel(i,j,k) < 200.)     &
        kntvel(k)=kntvel(k)+1
    ENDDO
    ENDDO
    ENDDO
    DEALLOCATE(hmin)
    DEALLOCATE(hmax)
    DEALLOCATE(itmp)
    print*,'finishing reading hdf'
    WRITE(6,'(a)') '  k   n ref    n vel'
    DO k=1,nz
      WRITE(6,'(i3,2i10)') k,kntref(k),kntvel(k)
    END DO
    istatus=1
  END IF
  RETURN
!
!  Destination for hdf read error
!

115   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in HDFREAD'

  istatus=-11
  RETURN

END SUBROUTINE read1rad
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GET_VIS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE get_vis(solar_alt,ni,nj,                                     & 1
           cloud_frac_vis_a,albedo,r_missing,istatus)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read satellite visible imagery data from LAPS ".lvd" file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!   04/1996   Based on the LAPS cloud analysis code of 1995.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!  INPUT:
  INTEGER :: ni,nj
  REAL :: solar_alt(ni,nj)
  REAL :: r_missing
!
!  OUTPUT:
!
  REAL :: albedo(ni,nj)
  REAL :: cloud_frac_vis_a(ni,nj)
!
!  LOCAL
!
  INTEGER :: ihist_alb(-10:20)
  INTEGER :: ihist_frac_sat(-10:20)
  INTEGER :: ih_cf_sat(-10:20)
!
!-----------------------------------------------------------------------
!
!  Misc. variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,iscr_alb
  INTEGER :: ibeg,iend,jbeg,jend
  INTEGER :: n_missing_albedo
  REAL :: albedo_to_cf,cloud_frac_vis
  REAL :: frac,term1,term2
  INTEGER :: istatus,iscr_frac_sat
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Initialize all histogram arrays.
!
!-----------------------------------------------------------------------
!
  DO i=-10,20
    ihist_alb(i) = 0
    ihist_frac_sat(i) = 0
    ih_cf_sat(i) = 0
  END DO
!
!
!-----------------------------------------------------------------------
!
!  Read the interpolated satellite albedo data.
!  These data were created by laps2arps.f (based on ter2arps)
!  by interpolating OLAPS (vortex95) grid satellite data onto
!  the ARPS grid.
!
!-----------------------------------------------------------------------
!
  ibeg=MAX(1,ni/2-5)
  iend=MIN(ni,ni/2+5)
  jbeg=MAX(1,nj/2-5)
  jend=MIN(nj,nj/2+5)
!
  WRITE(6,*) ' Satellite albedo value examples in the mid-domain:'

  WRITE(6,601) (i,i=ibeg,iend)
  601   FORMAT(/1X,3X,11(2X,i2,2X))
  WRITE(6,602) (j,(albedo(i,j),i=ibeg,iend)                             &
                ,j=jend,jbeg,-1)
  602   FORMAT(1X,i3,11F6.2)
!
  n_missing_albedo = 0
!
!-----------------------------------------------------------------------
!
!  Horizontal array loop
!
!-----------------------------------------------------------------------
!
  DO i = 1,ni
    DO j = 1,nj
!
!-----------------------------------------------------------------------
!
!  We now only use the VIS data if the solar alt exceeds 15 deg
!
!-----------------------------------------------------------------------
!
      IF (solar_alt(i,j) > 15.0 .AND. albedo(i,j) /= r_missing) THEN
!
!-----------------------------------------------------------------------
!
!  Translate the albedo into cloud fraction
!  Store histogram information for satellite data
!
!-----------------------------------------------------------------------
!
        iscr_alb  = nint(albedo(i,j)*10.)
        iscr_alb  = MIN(MAX(iscr_alb,-10),20)
        ihist_alb(iscr_alb) = ihist_alb(iscr_alb) + 1

        cloud_frac_vis = albedo_to_cf(albedo(i,j))
!
!-----------------------------------------------------------------------
!
!  Fudge the frac at low solar elevation angles
!  Note the ramp extrapolates down to 9 deg to account for slight
!  errors in determining the solar elevation
!
!-----------------------------------------------------------------------

        IF(solar_alt(i,j) < 20. .AND. solar_alt(i,j) >= 9.) THEN
          frac = (20. - solar_alt(i,j)) / 10.
          term1 = .13 * frac
          term2 = 1. + term1
          cloud_frac_vis = (cloud_frac_vis + term1) * term2
        END IF

        iscr_frac_sat = nint(cloud_frac_vis*10.)
        iscr_frac_sat = MIN(MAX(iscr_frac_sat,-10),20)
        ihist_frac_sat(iscr_frac_sat) =                                 &
                   ihist_frac_sat(iscr_frac_sat) + 1
!
!-----------------------------------------------------------------------
!
!  Make sure satellite cloud fraction is between 0 and 1
!
!-----------------------------------------------------------------------
!
        IF(cloud_frac_vis <= 0.0) cloud_frac_vis = 0.0
        IF(cloud_frac_vis >= 1.0) cloud_frac_vis = 1.0
        cloud_frac_vis_a(i,j) = cloud_frac_vis

        iscr_frac_sat = nint(cloud_frac_vis*10.)
        iscr_frac_sat = MIN(MAX(iscr_frac_sat,-10),20)
        ih_cf_sat(iscr_frac_sat) =                                      &
                   ih_cf_sat(iscr_frac_sat) + 1

      ELSE   ! albedo(i,j) = r_missing
        n_missing_albedo =  n_missing_albedo + 1
        cloud_frac_vis_a(i,j) = r_missing
      END IF  ! albedo(i,j).ne.r_missing

    END DO ! j
  END DO ! i

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

    IF(n_missing_albedo == ni*nj) THEN ! Return with status = 0
      WRITE(6,*)' All albedos were missing - return from get_vis'
      istatus = 0
      RETURN
    END IF

    WRITE(6,*)'       HISTOGRAMS'
    WRITE(6,*)'  I   Alb  CFsat  CFsato'
    DO i = -5,15
      WRITE(6,11)i,ihist_alb(i),ihist_frac_sat(i),ih_cf_sat(i)
      11      FORMAT(i4,i5,i7,i8)
    END DO ! i
!
    WRITE(6,*) ' ==== get_vis: Albedo derived cld cvr'
    WRITE(6,651) (i,i=ibeg,iend)
    651   FORMAT(/1X,3X,11(2X,i2,2X))
    WRITE(6,652) (j,(cloud_frac_vis_a(i,j),i=ibeg,iend)                 &
                  ,j=jend,jbeg,-1)
    652   FORMAT(1X,i3,11F6.2)

    WRITE(6,*)
    istatus = 1

    RETURN
END SUBROUTINE get_vis

!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION ALBEDO_TO_CF                   ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION albedo_to_cf( albedo )
!
!-----------------------------------------------------------------------
!
!  Old version assumes 0.15 land albedo and 0.85 cloud albedo if
!  the albedo fields are tuned properly.
!  cloud_frac_vis = (albedo - 0.15) / (0.85 - 0.15)
!
!  This version assumes 0.15 land albedo and 0.556 cloud albedo to
!  compensate for a bias in the albedo fields.
!  New value of "Cloud counts" should be 156 in normalize_laps.
!  to fix the problem.
!
!-----------------------------------------------------------------------
  REAL :: land_albedo,max_albedo
  PARAMETER (land_albedo=0.10,max_albedo=0.40)
  REAL :: rng_albedo
  PARAMETER (rng_albedo=(max_albedo-land_albedo))

  cloud_frac_vis = (albedo-land_albedo)/(max_albedo-land_albedo)

!-----------------------------------------------------------------------
!  Move out by 30% on either side
!-----------------------------------------------------------------------

!  cloud_frac_vis = 0.5 + (cloud_frac_vis - 0.5) * 1.7

!-----------------------------------------------------------------------
!  Add a 10% fudge factor to compensate for a drift in sensitivity
!-----------------------------------------------------------------------

!  cloud_frac_vis = cloud_frac_vis + 0.10

  albedo_to_cf = MIN(MAX(cloud_frac_vis,0.),1.0)

  RETURN
  END FUNCTION albedo_to_cf
!
!##################################################################
!##################################################################
!######                                                      ######
!######              SUBROUTINE SPREAD2                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE spread2 (cld_snd,wt_snd,i_snd,j_snd,n_cld_snd                & 10
           ,max_cld_snd,nz,iadas,jadas,k,cover,wt)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine inserts the cloud sounding into the analysis arrays
!  at one point location.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!  07/95
!
!  MODIFICATION HISTORY:
!  03/20/96. (Jian Zhang)
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: nz           ! number of cloud analysis levels
  INTEGER :: max_cld_snd  ! max. possible # of cloud soundings
                        ! in the domain

  INTEGER :: n_cld_snd    ! sequential index of cloud sounding sta.

  INTEGER :: iadas        ! i-loc of cld sounding stn in ADAS grid
  INTEGER :: jadas        ! j-loc of cld sounding stn in ADAS grid
  INTEGER :: k            ! the level where obs. cloud locates

  INTEGER :: i_snd (max_cld_snd)    ! array for i-location of
                                  ! cloud sounding stations
  INTEGER :: j_snd (max_cld_snd)    ! array for j-location of
                                  ! cloud sounding stations

  REAL :: cld_snd(max_cld_snd,nz)   ! cloud cover sounding array
  REAL :: wt_snd(max_cld_snd,nz)    ! weight array for cloud sounding

  REAL :: cover                     ! obs. cloud coverage
  REAL :: wt                        ! weight for the cloud cover obs.

  cld_snd(n_cld_snd,k) = cover
  wt_snd(n_cld_snd,k) = wt
  i_snd(n_cld_snd) = iadas
  j_snd(n_cld_snd) = jadas

  RETURN
END SUBROUTINE spread2
!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION T_GROUND_K                     ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION t_ground_k (t_sfc_k,solar_alt,solar_ha                       &
                ,solar_dec,rlat,cvr_snow,r_missing_data,i,j,ni,nj)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate the ground skin temperature as a function of
!  the surface air temperature and solar positions.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  (Jian Zhang)
!          Based on the LAPS cloud analysis code of 09/95
!
!  MODIFICATION HISTORY:
!  04/29/96  J. Zhang
!            Added ARPS format documents.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarition
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!  INPUT:
  INTEGER :: i,j,ni,nj
  REAL :: t_sfc_k            ! surface air temperature
  REAL :: solar_alt          ! solar altitude angle
  REAL :: solar_ha           ! solar hour angle
  REAL :: solar_dec          ! solar declination angle
  REAL :: rlat               ! latitude of the grid point
  REAL :: cvr_snow           ! ground snow cover at the grid point
  REAL :: r_missing_data     ! flag for data hole
!
!  OUTPUT:
  REAL :: t_ground_k         ! ground skin temperature
!
!  LOCAL:
  REAL :: high_alt,low_alt,corr_low,corr_high,corr_ramp
  REAL :: t_ground_fullcorr,t_snow_k,corr_negonly                       &
       ,solar_transit_alt,corr
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  solar_transit_alt = 90. - ABS(rlat - solar_dec)

  IF(solar_ha < 0.) THEN ! morning
    high_alt = solar_transit_alt * .66
    low_alt  = solar_transit_alt * .36
!       low_alt  = high_alt - 11.
  ELSE                     ! afternoon
    high_alt = solar_transit_alt * .79
    low_alt  = solar_transit_alt * .49
!       low_alt  = high_alt - 11.
  END IF

  corr_high = 4.
  corr_low = -4.
  corr_ramp = (corr_high - corr_low) / (high_alt - low_alt)
!
!-----------------------------------------------------------------------
!
!  Warmer at day, colder at night
!
!-----------------------------------------------------------------------
!
  IF(solar_alt < low_alt)THEN
    corr = corr_low

  ELSE IF(solar_alt > high_alt)THEN
    corr = corr_high

  ELSE ! ramp the function
    corr = corr_low + (solar_alt - low_alt) * corr_ramp

  END IF

  IF(i == ni/2 .AND. j == nj/2) THEN
    WRITE(6,600) i,j
    600     FORMAT(1X,'Example solar position parms for point (',       &
                      2I2,')')
    WRITE(6,601)
    601     FORMAT(1X,'  hangle  altitude transalt  ratio  ',           &
                   ' highalt   lowalt  correctn')
    WRITE(6,602) solar_ha,solar_alt,solar_transit_alt                   &
               ,solar_alt/solar_transit_alt                             &
               ,high_alt,low_alt,corr
    602     FORMAT(1X,7F9.2)
  END IF

  corr_negonly = MIN(corr,0.0)            ! Only colder at night

  IF(cvr_snow /= r_missing_data) THEN
    t_ground_fullcorr = t_sfc_k + corr  ! Warmer at day, colder at night
    t_snow_k = t_sfc_k + corr_negonly   ! Only colder at night
    t_snow_k = MIN(t_snow_k,273.15)     ! Snow can't be above 0C
    t_ground_k = t_ground_fullcorr * (1.-cvr_snow)                      &
                          + t_snow_k * cvr_snow
  ELSE
    t_ground_k = t_sfc_k + corr_negonly ! Only colder at night
  END IF

  RETURN
  END FUNCTION t_ground_k
!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION VASRAD                         ######
!######                                                      ######
!##################################################################
!##################################################################
!


  REAL FUNCTION vasrad(tau,temp,tsfc,kchan,nl,emiss)
!dis
!dis
!dis
!dis    Forecast Systems Laboratory
!dis    NOAA/OAR/ERL/FSL
!dis    325 Broadway
!dis    Boulder, CO     80303
!dis
!dis    Forecast Research Division
!dis    Local Analysis and Prediction Branch
!dis    LAPS
!dis
!dis    This software and its documentation are in the public domain and
!dis    are furnished "as is."  The United States government, its
!dis    instrumentalities, officers, employees, and agents make no
!dis    warranty, express or implied, as to the usefulness of the software
!dis    and documentation for any purpose.  They assume no responsibility
!dis    (1) for the use of the software and documentation; or (2) to provide
!dis     technical support to users.
!dis
!dis    Permission to use, copy, modify, and distribute this software is
!dis    hereby granted, provided that the entire disclaimer notice appears
!dis    in all copies.  All modifications to this software must be clearly
!dis    documented, and are solely the responsibility of the agent making
!dis    the modifications.  If significant modifications or enhancements
!dis    are made to this software, the FSL Software Policy Manager
!dis    (softwaremgr@fsl.noaa.gov) should be notified.
!dis
!dis
!dis
!dis
!    routine received 5 aug 1993 (from wisconsin) as a replacement for
!    an earlier version that was out of date. and sent by mistake

! $ function vasrad(tau,temp,tsfc,kchan,nl,emiss)   (btr)
! $ vasrad - get vas radiance from temperature profile
! $ input:
! $     tau = (r) array of vas transmittances
! $     temp = (r) array of atmospheric temperatures
! $     tsfc = (r) temperature of surface
! $     kchan = (i) channel number
! $     nl = (i) level number of surface
! $     emiss = (r) surface emissivity
! $ output description:
! $     vas radiance
! $$ vasrad = sounder, vas, convert
  COMMON/gde/gv(12),dv(12),ev(12)
  DIMENSION tau(*),temp(*)
  t1=temp(1)
  b1=vplanc(t1,kchan)
  tau1=tau(1)
  rad=0.
  DO i=2,nl
    t2=temp(i)
    b2=vplanc(t2,kchan)
    tau2=tau(i)
    rad=rad+.5*(b1+b2)*(tau1-tau2)
    b1=b2
    tau1=tau2
  END DO
  bs=vplanc(tsfc,kchan)
  rad=rad+emiss*bs*tau(nl)
  rad=rad+dv(kchan)
  rad=AMAX1(rad,.001)
  vasrad=rad
  RETURN
  END FUNCTION vasrad

!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION VPLANC                         ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION vplanc(t,k)

  COMMON/plankv/fnu(12),fk1(12),fk2(12),tc(2,12)
  tt=tc(1,k)+tc(2,k)*t
  expn=EXP(fk2(k)/tt) - 1.
  vplanc=fk1(k)/expn
  RETURN
  END FUNCTION vplanc
!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION VBRITE                         ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION vbrite(r,k)

  COMMON/plankv/fnu(12),fk1(12),fk2(12),tc(2,12)
  expn=fk1(k)/r+1.
  tt=fk2(k)/ALOG(expn)
  vbrite=(tt-tc(1,k))/tc(2,k)
  RETURN
  END FUNCTION vbrite

!
!##################################################################
!##################################################################
!######                                                      ######
!######              FUNCTION VDBDTB                         ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION vdbdtb(tbb,k)

  COMMON/plankv/fnu(12),fk1(12),fk2(12),tc(2,12)
  t=tbb
  tt=tc(1,k)+tc(2,k)*t
  expn=EXP(fk2(k)/tt) - 1.
  b=fk1(k)/expn
  bb=b*b
  tt=t*t
  expn=EXP(fk2(k)/t)
  q=fk2(k)/fk1(k)
  vdbdtb=q*expn*bb/tt
  RETURN
  END FUNCTION vdbdtb


SUBROUTINE vasrte(tau,temp,tsfc,rad,dbdt,tbb,dbdtbb,kchan,lgnd)
  COMMON/plankv/fnu(12),fk1(12),fk2(12),tc(2,12)
  COMMON/gde/gv(12),dv(12),ev(12)
  DIMENSION tau(*),temp(*),dbdt(*)
  DATA nl/40/
  f1=fk1(kchan)
  f2=fk2(kchan)
  q=f2/f1
  t1=temp(1)
  b1=vplanc(t1,kchan)
  bb=b1*b1
  tt=t1*t1
  ex=q*EXP(f2/t1)
  dbdt(1)=ex*bb/tt
  tau1=tau(1)
  rad=0.
  DO i=2,nl
    t2=temp(i)
    b2=vplanc(t2,kchan)
    bb=b2*b2
    tt=t2*t2
    ex=q*EXP(f2/t2)
    dbdt(i)=ex*bb/tt
    tau2=tau(i)
    IF(i > lgnd) GO TO 100
    rad=rad+.5*(b1+b2)*(tau1-tau2)
!
!            02-may-84       debug
!
!            type 1000, i, b1, b2, tau1, tau2, rad
!    1000            FORMAT (i3, 5(1X, e12.5))
    100 b1=b2
    tau1=tau2
  END DO
  bs=vplanc(tsfc,kchan)
  rad=rad+bs*tau(lgnd)
  rad=rad+dv(kchan)
  rad=AMAX1(rad,.001)
  tbb=vbrite(rad,kchan)
  tt=tbb*tbb
  bb=rad*rad
  ex=q*EXP(f2/tbb)
  dbdtbb=ex*bb/tt
  RETURN
END SUBROUTINE vasrte


SUBROUTINE plnkiv(isat,calib_fname) 2,3
!
!         1 7-apr-84 change in file name specification
!  real*8 file
!
  COMMON/plankv/fnu(12),fk1(12),fk2(12),tc(2,12)
  COMMON/gde/gv(12),dv(12),ev(12)
  COMMON/tsurfc/tsc(4,2)
  COMMON/tsurfe/tse(4,2)
  COMMON/use/iuch(12,2)
  COMMON/chkcof/chkc(5)
  DIMENSION buf(506)
!
!         1 7-apr-84 change in file name specification
!  data file/'vastaucf'/,len/506/
!
  CHARACTER (LEN=132) :: calib_fname
!   byte file(28)
!   data file/'',
!  1   'v','a','s','t','a','u','c','f','.','d','a','t',"0/
  DATA LEN/506/
  DATA iu/13/,nc/12/
!
!         1 3-jun-j4 debug
!
!            type 1090, isat
!1090            FORMAT (/, 1X, ' satellite ', i3)
  CALL dopen(calib_fname,iu,LEN)
  irec=isat*(nc+1)
  CALL dread(iu,irec,buf)
  CALL dclose(iu)
  mc=nc*2
  DO j=1,nc
    fnu(j)=buf(j)
    fk1(j)=buf(j+nc)
    fk2(j)=buf(j+mc)
  END DO
!
!         1 3-jun-84 debug
!
!            type 1000, (fnu(iii), fk1(iii), fk2(iii), iii = 1, nc)
!  1000            FORMAT (/, 1X, 'fnu         fk1          fk2', /,     &
!                        (3E12.4))
  n=nc*3
  DO j=1,nc
    DO i=1,2
      n=n+1
      tc(i,j)=buf(n)
    END DO
  END DO
!
!         1 3-jun-84 debug
!
!            type 1010, ((tc(iii, jjj), iii = 1, 2), jjj = 1, nc)
!  1010            FORMAT (/, 1X, 'tc1        tc2', /, (2E12.4))
  DO i=1,nc
    n=n+1
    gv(i)=buf(n)
  END DO
!
!         1 3-jun-84 debug
!
!            type 1020, (gv(iii), iii = 1, nc)
!  1020            FORMAT (/, 1X, 'gb', /, (e12.4))
  DO i=1,nc
    n=n+1
    dv(i)=buf(n)
  END DO
!
!         1 3-jun-84 debug
!
!            type 1030, (dv(iii), iii = 1, nc)
!  1030            FORMAT (/, 1X, 'dv', /, (e12.4))
  DO i=1,nc
    n=n+1
    ev(i)=buf(n)
  END DO
!
!         1 3-jun-84 debug
!
!            type 1040, (ev(iii), iii = 1, nc)
!  1040            FORMAT (/, 1X, 'ev', /, (e12.4))
  n=100
  DO j=1,2
    DO i=1,4
      n=n+1
      tsc(i,j)=buf(n)
    END DO
  END DO
!
!         1 3-jun-84 debug
!
!            type 1050, ((tsc(iii, jjj), iii = 1, 2), jjj = 1, 4)
!  1050            FORMAT (/, 1X, 'tsc1          tsc2', /,               &
!                        (2E12.4))
  n=110
  DO j=1,2
    DO i=1,nc
      n=n+1
      iuch(i,j)=buf(n)
    END DO
  END DO
!
!         1 3-jun-84 debug
!
!            type 1060, ((iuch(iii, jjj), jjj = 1, 2), iii = 1, nc)
!  1060            FORMAT (/, 1X, 'iuch1         iuch2', /,              &
!                        (2I12))
  n=140
  DO i=1,5
    n=n+1
    chkc(i)=buf(n)
  END DO
  n=150
!
!         1 3-jun-84 deubg
!
!            type 1070, (chkc(iii), iii = 1, 5)
!  1070            FORMAT (/, 1X, 'chkc', /, (5E12.4))
  DO j=1,2
    DO i=1,4
      n=n+1
      tse(i,j)=buf(n)
    END DO
  END DO
!
!         1 3-jun-84 debug
!
!            type 1080, ((tse(iii, jjj), jjj = 1, 2), iii = 1, 4)
!  1080            FORMAT (/, 1X, 'tse1         tse2', /,                &
!                        (2E12.4))
  RETURN
END SUBROUTINE plnkiv


SUBROUTINE dopen (FILE, iu, LEN) 1
!
!         1 7-apr-84 subroutine added
!
  COMMON/dio/lenn
!    byte file(1)
  CHARACTER (LEN=132) :: FILE
!
  PRINT *, 'Opening calibration file=',FILE
  lenn = LEN
  OPEN (UNIT=iu, FILE=trim(FILE),                                       &
      FORM='formatted')
!
  RETURN
END SUBROUTINE dopen


SUBROUTINE dread (iu, irec, buf) 1
!
!         1 7-apr-84 subroutine added
!
!wdt update
  IMPLICIT NONE
  INTEGER iu,irec
  INTEGER lenn
  COMMON/dio/lenn
  REAL buf(lenn), allbuf(506,39)
  INTEGER :: test
  DATA test /0/
  !wdt update
  SAVE allbuf
  INTEGER :: i,j
!
!    read (unit=iu, rec=irec) (buf(i), i = 1, lenn)
  IF (test == 1) THEN
    CONTINUE
  ELSE
    test = 1
    PRINT *, 'Reading in calibration file.'
    READ (iu, *) ((allbuf(i,j), i=1,506),j=1,39)
    PRINT *, 'Done reading calibration file.'
  END IF

!    go here directly if the file has been read before
  DO i = 1,lenn
    buf(i) = allbuf(i,irec)
  END DO

!
  RETURN
END SUBROUTINE dread


SUBROUTINE dclose (iu) 1
!
!         1 7-apr-84 subroutine added
!
  COMMON/dio/lenn
  CLOSE (UNIT=iu)
  RETURN
END SUBROUTINE dclose