subroutine reorg_metar_cloud(cdata,nreal,ndata,cdata_all,maxobs,ngrid)
!$$$ subprogram documentation block
! . . . .
! subprogram: reorg_metar_cloud define a closest METAR cloud observation for each grid point
!
! PRGMMR: Ming Hu ORG: GSD/AMB DATE: 2009-09-21
!
! ABSTRACT:
!
! PROGRAM HISTORY LOG:
! 2010-04-21 Hu initial
! 2011-10-21 Hu add code to delete dust and haze report
! 2015-08-06 S.Liu cross-check low cloud and T-Td and reject
! low cloud if T-Td > 2.0 deg
! 2016-06-21 S.Liu give number precision
!
! input argument list:
! cdata - METRA cloud observation
! nreal - first dimension of cdata
! ndata - second dimension of cdata
! maxobs - maximum number of cdata_all
! ndata - number of type "obstype" observations retained for further processing
!
! output argument list:
! cdata_all - reorganized METAR observation
! ngrid - number of type metar cloud observations reorganized to the grid
!
! USAGE:
! INPUT FILES:
! OUTPUT FILES:
!
! REMARKS:
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90
! MACHINE: Linux cluster (WJET)
!
!$$$
!
!_____________________________________________________________________
!
use kinds, only: r_kind,i_kind,r_double
use gridmod, only: nlon,nlat
use constants, only: one,half
use rapidrefresh_cldsurf_mod, only: metar_impact_radius
implicit none
! Declare passed variables
integer(i_kind) ,intent(in) :: nreal
integer(i_kind) ,intent(in) :: ndata
integer(i_kind) ,intent(in) :: maxobs
real(r_kind),dimension(nreal,ndata) ,intent(inout) :: cdata
real(r_kind),dimension(nreal,maxobs) ,intent(out):: cdata_all
integer(i_kind) ,intent(out):: ngrid
! local variable
!
integer(i_kind) :: ista_prev,ista_prev2,ista_save
INTEGER(i_kind),allocatable :: first_sta(:,:)
INTEGER(i_kind),allocatable :: next_sta(:)
INTEGER(i_kind) :: null_p
PARAMETER ( null_p = -1 )
INTEGER(i_kind) :: nsta_cld
INTEGER(i_kind),allocatable :: sta_cld(:)
integer(i_kind) :: isprd,isprd2,iout
integer(i_kind) :: aninc_cld_p
integer(i_kind) :: i,j,k,i1,j1,ic,ista
integer(i_kind) :: ista_min
real(r_kind) :: min_dist, dist
real(r_kind) :: awx, cg_dewpdep, cldamt,cldhgt,cl_base_ista
real(r_kind) :: cl_mt_ista
integer(i_kind) :: idust,ihaze
real(r_double) rstation_id
character(8) :: cstation1,cc,ci
equivalence(cstation1,rstation_id)
real(r_kind) :: spval_p
parameter (spval_p = 99999._r_kind)
!
!
isprd=int(metar_impact_radius + half)
if(isprd <= 0) return
if(ndata <= 0) return
ngrid = 0
aninc_cld_p = isprd
isprd2 = isprd * isprd
!
! check duplicate observation and pick the one near analyis time
!
DO i=1,ndata-1
DO ic = i+1, ndata
rstation_id=cdata(1,ic)
cc=cstation1
rstation_id=cdata(1,i)
ci=cstation1
if(trim(cc) == trim(ci) ) then
if(abs(cdata(21,ic)) >= abs(cdata(21,i))) then ! 21= observation time
cdata(22,ic)=100.0_r_kind ! 22= usage
else
cdata(22,i)=100.0_r_kind
endif
endif
ENDDO ! ic
ENDDO ! i
!
! Check Haze and Dust station
!
DO i=1,ndata
rstation_id=cdata(1,i)
cc=cstation1
if(cdata(22,i) < 50 ) then
idust=0
ihaze=0
DO j=1,3
awx = cdata(17+j,i) ! weather
if(awx==5._r_kind .or. awx==105._r_kind .or. awx==104._r_kind) ihaze=ihaze+1
if( (awx >=6._r_kind .and. awx <= 9._r_kind) .or. &
(awx >=30._r_kind .and. awx <= 35._r_kind) .or. &
(awx ==76._r_kind .or. awx == 111._r_kind) &
) idust=idust+1
ENDDO ! j
!* check cloud base and T-Td
cl_base_ista=spval_p
DO j=1,3
cldamt = cdata(5+j,i) ! cloud amount
cldhgt = int(cdata(11+j,i)) ! cloud bottom height
if(abs(cldamt) < spval_p .and. abs(cldhgt) < spval_p) then
if(cl_base_ista>cldhgt)then
if( (cldamt > 3.5_r_kind .and. cldamt < 9.5_r_kind ) .or. abs(cldamt-12._r_kind) < 0.0001_r_kind ) then
cl_base_ista = cldhgt
cl_mt_ista = cldamt
endif ! cloud ceiling
endif
endif
ENDDO ! j
!* case1: low cloud with full coverage, but saturation height derived from T-Td(>1.5 deg) > 150
if(cl_base_ista < 100.0_r_kind .and. &
cl_mt_ista>6.0_r_kind .and. cl_mt_ista < 9.5_r_kind .and. &
cdata(24,i)>2.5 ) then
cdata(22,i)=200
write(6,*)'block station::',cc
end if
if(ihaze > 0 .or. idust > 0 ) then
! write(*,'(a,a8,3F10.1)') 'dust or haze report at station ',rstation_id, (cdata(17+j,i),j=1,1)
cg_dewpdep = 100._r_kind * cdata(24,i) ! 100. = dry adiabatic lapse rate in m/K
DO j=1,3
cldamt = cdata(5+j,i) ! cloud amount
cldhgt = int(cdata(11+j,i)) ! cloud bottom height
if(abs(cldamt) < spval_p .and. abs(cldhgt) < spval_p) then
if( (cldamt > 3.5_r_kind .and. cldamt < 9.5_r_kind ) .or. abs(cldamt-12._r_kind) < 0.0001_r_kind ) then
cl_base_ista = cldhgt
! Check - compare estimated cloud base from dewpoint depression (calc above as cg_dewpdep)
! with ceiling height AGL. We allow a slop room of 200m extra. Stan B. and John B.-31Aug2011
if (cg_dewpdep > -10._r_kind .and. (cg_dewpdep-300._r_kind > cl_base_ista )) then
cdata(22,i)=200
write (6,'(a,a,3(a,f10.1))') 'Dust-ceiling IDed ', rstation_id, &
' dewpoint depr *100=',cg_dewpdep, 'ceiling=',cl_base_ista
end if
if (cg_dewpdep < -990._r_kind) then
cdata(22,i)=200
write (6,*) 'Dust-ceiling possible, could not calc dewpoint, do not use ceiling ',rstation_id
end if
endif ! cloud ceiling
endif
ENDDO ! j
endif ! ihaze > 0 .or. idust > 0
endif
ENDDO ! i
!
!
! find the nearest METAR station for each grid point
!
allocate(first_sta(nlon,nlat))
allocate(next_sta(ndata))
allocate(sta_cld(ndata))
first_sta= null_p
next_sta = null_p
sta_cld = 0
! -- Then fill arrays based on where obs are
! -- NEXT_STA - pointer to next ob found in grid volume
! from the previous ob in that volume
!-----------------------------------------------------------
do ista = 1,ndata
if (cdata(22,ista) < 50) then
if (cdata(6,ista) >= -one .and. cdata(6,ista) < 100.0_r_kind) then
! when cloud data are available
i = int(cdata(2,ista))
j = int(cdata(3,ista))
if (first_sta(i,j) == null_p) then
first_sta(i,j) = ista
else
ista_prev = first_sta(i,j)
do while (ista_prev /= null_p )
ista_prev2= next_sta(ista_prev)
ista_save = ista_prev
ista_prev = ista_prev2
enddo
next_sta(ista_save) = ista
end if
endif
endif
enddo
iout=0
DO j = 1,nlat,aninc_cld_p
DO i = 1,nlon,aninc_cld_p
nsta_cld = 0
!sb -- Find all stations w/ cloud data within encompassing box
!mh The box is decide by the influence radius of the analysis
do j1 = max(1,j-isprd),min(nlat,j+aninc_cld_p+isprd)
do i1 = max(1,i-isprd),min(nlon,i+aninc_cld_p+isprd)
if (first_sta(i1,j1) /= null_p) then
ista_prev = first_sta(i1,j1)
do while (ista_prev /= null_p )
nsta_cld = nsta_cld + 1
sta_cld(nsta_cld) = ista_prev
ista_prev2= next_sta(ista_prev)
ista_prev = ista_prev2
enddo
end if
enddo ! j1
enddo ! i1
!
!sb - This is the big grid pt loop. Walk thru each
! individual grid points within 10x10 box.
if(nsta_cld > 0 ) then
do j1 = j,min(nlat,j+aninc_cld_p - 1)
do i1 = i,min(nlon,i+aninc_cld_p - 1)
!sb - Find closest cloud station to grid point
min_dist = 1.e10_r_kind
do ic= 1,nsta_cld
ista = sta_cld(ic)
dist = (float(i1)-cdata(2,ista))*(float(i1)-cdata(2,ista)) &
+(float(j1)-cdata(3,ista))*(float(j1)-cdata(3,ista))
if (dist < min_dist .and. dist < float(isprd2)) then
min_dist = dist
ista_min = ista
end if
end do ! ic find the closest cloud station
if (min_dist < 1.e9_r_kind) then
if (i1 > 1 .and. i1 < nlon .and. j1 > 1 .and. j1 < nlat) then
iout = iout + 1
if(iout > maxobs) then
write(6,*)'reorg_metar_cloud: ***Error*** ndata > maxobs '
call stop2(50)
end if
do k=1,nreal
cdata_all(k,iout) = cdata(k,ista_min)
enddo
cdata_all(24,iout) = cdata_all(2,iout) ! save observaion station i
cdata_all(25,iout) = cdata_all(3,iout) ! save observaion station j
cdata_all(2,iout) = float(i1) ! grid index i
cdata_all(3,iout) = float(j1) ! grid index j
cdata_all(23,iout)=sqrt(min_dist) ! distance from station
endif
endif
enddo ! j1
enddo ! i1
endif ! nsta_cld > 0
enddo ! j
enddo ! i
ngrid = iout
deallocate(first_sta)
deallocate(sta_cld)
deallocate(next_sta)
end subroutine reorg_metar_cloud