subroutine map_ctp_lar(ib,jb,nx,ny,nn_obs,numsao,data_s,sat_ctp,sat_tem,w_frac,w_lwp,nlev_cld)
!
!$$$ subprogram documentation block
! . . . .
! subprogram: map_ctp_lar map Langley cloud product to analysis grid
!
! PRGMMR: Shun Liu ORG: GSD/AMB DATE: 2006-03_10
!
! ABSTRACT:
! This subroutine map Langley cloud product to analysis grid, copy from map_ctp
!
! PROGRAM HISTORY LOG:
! 2009-01-20 Hu Add NCO document block
!
!
! input argument list:
! ib - begin i point of this domain
! jb - begin j point of this domain
! nx - no. of lons on subdomain (buffer points on ends)
! ny - no. of lats on subdomain (buffer points on ends)
! nn_obs - 1st dimension of observation arry data_s
! numsao - number of observation
! data_s - observation array for GOES cloud products
!
! output argument list:
! sat_ctp - GOES cloud top pressure in analysis grid
! sat_tem - GOES cloud top temperature in analysis grid
! w_frac - GOES cloud coverage in analysis grid
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
!
! adapted according to RUC subroutine rd_cld
! *
! * This routine reads NESDIS (Madison, WI) cloud product produced
! * from GOES sounder data. The original product is reprocessed onto
! * MAPS40 grid boxes. There could be more than one cloud product
! * in a grid-box, so we use the nearest one that falls in the
! * grid. The routine combines GOES-8 and 10 products.
!
! ===== History =====
!
! * Internal variables:
! CTP_E, CTP_W Soft-linked filename for ascii GOES Clouds
!
! * Working variables:
!
! * Working variables used for sorting max size of 10:
! Pxx, Txx, xdist,xxxdist (R4)
! Fxx, Nxx, index, jndex (I4)
! ioption (I4) = 1 if selection is nearest neighbor
! = 2 if selection is median of samples
!
!
! * Output variables on gridpoint (Nx,Ny):
! sat_ctp, sat_tem (R4) Cloud-top pressure and temperature
! w_frac (R4) Effective fractional cloud coverage, option=1
! fractional coverage within RUC grid, option=2
! w_eca (R4) Effective fractional cloud regardless option
! (effective cloud amount - eca)
! nlev_cld (I4) Number of cloud levels. TO BE USED LATER
! to incorporate multi-level cloud
!
! * Calling routines
! sorting
! sortmed
!
! *
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90
! MACHINE: Linux cluster (WJET)
!
!$$$
!
!_____________________________________________________________________
!
use kinds, only: r_kind,r_single,i_kind
use constants, only: zero,one_tenth,one,deg2rad
implicit none
! input-file variables:
INTEGER(i_kind),intent(in) :: Nx, Ny
INTEGER(i_kind),intent(in) :: ib, jb
INTEGER(i_kind),intent(in) :: numsao, nn_obs
real(r_kind),dimension(nn_obs,numsao):: data_s
! Output
real(r_single), intent(out) :: sat_ctp(Nx,Ny)
real(r_single), intent(out) :: sat_tem(Nx,Ny)
real(r_single), intent(out) :: w_lwp(Nx,Ny)
real(r_single), intent(out) :: w_frac(Nx,Ny)
!
! misc
integer(i_kind) :: nfov
parameter (nfov=650)
character header*80
! Working
real(r_kind) :: Pxx(Nx,Ny,nfov),Txx(Nx,Ny,nfov)
real(r_kind) :: PHxx(Nx,Ny,nfov),WPxx(Nx,Ny,nfov)
real(r_kind) :: xdist(Nx,Ny,nfov), xxxdist(nfov)
real(r_kind) :: fr,sqrt, qc, type
integer(i_kind) :: Nxx(Nx,Ny,nfov),index(Nx,Ny), jndex(nfov)
integer(i_kind) :: ioption
integer(i_kind) :: ipt,ixx,ii,jj,i,med_pt,igrid,jgrid &
,ncount,ncount1,ncount2,ii1,jj1,nobs,n
real(r_kind) :: xc
real(r_kind) :: yc
real(r_single) :: w_eca(Nx,Ny)
integer(i_kind) :: nlev_cld(Nx,Ny)
integer(i_kind) :: ios,cfov
!
! * Initialize outputs since GOES sounder do not scan all MAPS domain
!
do jj=1,Ny
do ii=1,Nx
sat_ctp (ii,jj) =-99999._r_kind
sat_tem (ii,jj) =-99999._r_kind
w_lwp (ii,jj) =-99999._r_kind
w_frac (ii,jj) =-99999._r_kind
nlev_cld (ii,jj) =-99999
index(ii,jj) = 0
enddo
enddo
! -- set ios as failed unless valid data points are found below
ios = 0
! -----------------------------------------------------------
! -----------------------------------------------------------
! Map each FOV onto RR grid points
! -----------------------------------------------------------
! -----------------------------------------------------------
do ipt=1,numsao
xc=data_s(2,ipt) - ib + 1.0_r_kind
yc=data_s(3,ipt) - jb + 1.0_r_kind
! write(6,*)'sat_tem::',data_s(2,ipt),data_s(3,ipt),ib,jb
if(data_s(8,ipt) > 650 ) cycle
! * XC,YC should be within subdomain boundary, i.e., XC,YC >0
if(XC >= 1._r_kind .and. XC < Nx .and. &
YC >= 1._r_kind .and. YC < Ny) then
ii1 = int(xc+0.5_r_kind)
jj1 = int(yc+0.5_r_kind)
do jj = max(1,jj1-2), min(ny,jj1+2)
if (jj1-1 >= 1 .and. jj1+1 <= ny) then
do ii = max(1,ii1-2), min(nx,ii1+2)
if (ii1-1 >= 1 .and. ii1+1 <= nx) then
! * We check multiple data within gridbox
if (index(ii,jj) < nfov) then
index(ii,jj) = index(ii,jj) + 1
Pxx(ii,jj,index(ii,jj)) = data_s(4,ipt)
Txx(ii,jj,index(ii,jj)) = data_s(5,ipt)
PHxx(ii,jj,index(ii,jj)) = data_s(6,ipt)
WPxx(ii,jj,index(ii,jj)) = data_s(7,ipt)
!mhu Nxx(ii,jj,index(ii,jj)) = int(data_s(5,ipt))
!mhu no cloud amount available, assign to 100
! Nxx(ii,jj,index(ii,jj)) = 100
nlev_cld(ii,jj) = 1
! write(6,*)'sat_tem1::',index(ii,jj),data_s(4,ipt),data_s(5,ipt),data_s(6,ipt),data_s(7,ipt)
xdist(ii,jj,index(ii,jj)) = sqrt( &
(XC+1-ii)**2 + (YC+1-jj)**2)
end if
endif
enddo ! ii
endif
enddo ! jj
endif ! observation is in the domain
enddo ! ipt
!
! * ioption = 1 is nearest neighrhood
! * ioption = 2 is median of cloudy fov
ioption = 2
!
do jj = 1,Ny
do ii = 1,Nx
if (index(ii,jj) < 3 ) then
! sat_ctp(ii,jj) = Pxx(ii,jj,1)
! sat_tem(ii,jj) = Txx(ii,jj,1)
! w_frac(ii,jj) = float(Nxx(ii,jj,1))/100.
! w_eca(ii,jj) = float(Nxx(ii,jj,1))/100.
elseif(index(ii,jj) >= 3) then
! * We decided to use nearest neighborhood for ECA values,
! * a kind of convective signal from GOES platform...
do i=1,index(ii,jj)
jndex(i) = i
xxxdist(i) = xdist(ii,jj,i)
enddo
call sorting(xxxdist,index(ii,jj),jndex)
! w_eca(ii,jj) = float(Nxx(ii,jj,jndex(1)))/100._r_kind
! * Sort to find closest distance if more than one sample
! if(ioption == 1) then !nearest neighborhood
! do i=1,index(ii,jj)
! jndex(i) = i
! xxxdist(i) = xdist(ii,jj,i)
! enddo
! call sorting(xxxdist,index(ii,jj),jndex)
! sat_ctp(ii,jj) = Pxx(ii,jj,jndex(1))
! sat_tem(ii,jj) = Txx(ii,jj,jndex(1))
! w_frac(ii,jj) = float(Nxx(ii,jj,jndex(1)))/100._r_kind
! endif
! * Sort to find median value
if(ioption == 2) then !pick median
do i=1,index(ii,jj)
jndex(i) = i
xxxdist(i) = Pxx(ii,jj,i)
enddo
call sortmed(xxxdist,index(ii,jj),jndex,fr)
med_pt = index(ii,jj)/2 + 1
sat_ctp(ii,jj) = Pxx(ii,jj,jndex(med_pt))
sat_tem(ii,jj) = Txx(ii,jj,jndex(med_pt))
w_lwp(ii,jj) = WPxx(ii,jj,jndex(med_pt))
if (sat_ctp(ii,jj).eq.-20) then
sat_ctp(ii,jj) = 1013. ! hPa - no cloud
w_frac(ii,jj)=0.0
nlev_cld(ii,jj) = 0
end if
!
! cloud fraction based on phase (0 are clear), what about -9 ????
if( sat_ctp(ii,jj) < 1012.99) then
cfov = 0
do i=1,index(ii,jj)
if(PHxx(ii,jj,i) .gt. 0.1) cfov = cfov + 1
enddo
w_frac(ii,jj) = float(cfov)/(max(1,index(ii,jj))) ! fraction
if( w_frac(ii,jj) > 0.01 ) nlev_cld(ii,jj) = 1
endif
! write(6,*)'sat_tem2::',index(ii,jj),sat_ctp(ii,jj),sat_tem(ii,jj)
endif
endif
enddo !ii
enddo !jj
return
end subroutine map_ctp_lar