subroutine read_NASA_LaRC_cloud(nread,ndata,nouse,infile,obstype,lunout,sis,nobs)
!
! PRGMMR: Shun Liu ORG: EMC DATE: 2013-05-14
!
! ABSTRACT:
! This routine read in NASA LaRC cloud products based on GSD's cloud processing
! code
!
! PROGRAM HISTORY LOG:
! 2014-12-03 derber - remove unused variables
!
! variable list
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90 + EXTENSIONS
! MACHINE: CCS
!
!$$$
!
!_____________________________________________________________________
!
use kinds, only: r_kind,i_kind,r_single
use constants, only: zero,deg2rad
use gridmod, only: regional,nlat,nlon,tll2xy,rlats,rlons
use mpimod, only: npe
implicit none
! Declare passed variables
character(len=*),intent(in ) :: obstype
character(len=20),intent(in ) :: sis
character(15), intent(in ) :: infile
integer(i_kind) ,intent(in ) :: lunout
integer(i_kind) ,intent(inout) :: nread,ndata,nouse
integer(i_kind) ,dimension(npe),intent(inout) :: nobs
! real(r_kind),dimension(nlat,nlon,nsig),intent(in):: hgtl_full
! Declare local parameters
integer(i_kind),parameter:: maxdat=8
! Declare local variables
logical outside !,good0,lexist1,lexist2,lexist3
real(r_kind),allocatable,dimension(:,:):: cdata_all
real(r_kind) dlat_earth,dlon_earth
real(r_kind) dlat,dlon
real(r_kind) dlatmax,dlonmax,dlatmin,dlonmin
real(r_kind) usage
integer(i_kind) nreal,nchanl,ilat,ilon
real(r_kind),parameter:: r360=360.0_r_kind
CHARACTER*80 satfile
! ****VARIABLES FOR THIS NETCDF FILE****
!
REAL(r_single), allocatable :: lat_l(:)
REAL(r_single), allocatable :: lon_l(:)
REAL(r_single), allocatable :: lwp_l(:)
REAL(r_single), allocatable :: teff_l(:)
REAL(r_single), allocatable :: ptop_l(:)
integer(i_kind), allocatable :: phase_l(:)
!
! array for RR
!
integer :: nfov
parameter (nfov=650)
integer i,k
integer(i_kind) :: east_time, west_time
integer :: maxobs,numobs
character*10 atime
!**********************************************************************
!
! END OF DECLARATIONS....start of program
!
! set geogrid fle name
!
! allocate (Pxx(nlon,nlat,nfov),Txx(nlon,nlat,nfov),WPxx(nlon,nlat,nfov))
! allocate (xdist(nlon,nlat,nfov), xxxdist(nfov))
! allocate (PHxx(nlon,nlat,nfov),index(nlon,nlat), jndex(nfov))
! index=0
!
! read in the NASA LaRC cloud data
write(6,*)'start to read NASA_LaRC_cloud::'
!maxobs=(1800*700 + 1500*850)*1
maxobs=(12250000)
allocate(cdata_all(maxdat,maxobs))
!satfile='lgycldbufr'
satfile=trim(infile)
allocate(lat_l(maxobs))
allocate(lon_l(maxobs))
allocate(ptop_l(maxobs))
allocate(teff_l(maxobs))
allocate(phase_l(maxobs))
allocate(lwp_l(maxobs))
usage=0
ptop_l=-9.0_r_kind
teff_l=-9.0_r_kind
lat_l =-9.0_r_kind
lon_l =-9.0_r_kind
lwp_l =-9.0_r_kind
phase_l=-9
call read_NASALaRC_cloud_bufr_survey(satfile,east_time, west_time)
call read_NASALaRC_cloud_bufr(satfile,atime, &
maxobs,numobs, ptop_l, teff_l, phase_l, lwp_l,lat_l, lon_l)
! write(6,*)'LaRC ptop =', (ptop_l(j),j=1,numobs,5000)
! write(6,*)'LaRC teff =', (teff_l(j),j=1,numobs,5000)
! write(6,*)'LaRC lat =', (lat_l(j),j=1,numobs,5000)
! write(6,*)'LaRC lon =', (lon_l(j),j=1,numobs,5000)
! write(6,*)'LaRC lwp =', (lwp_l(j),j=1,numobs,5000)
! write(6,*)'LaRC phase =', (phase_l(j),j=1,numobs,5000)
do i=1,numobs
dlat_earth=lat_l(i) !station lat (degrees)
dlon_earth=lon_l(i) !station lon (degrees)
if (dlon_earth>=r360) dlon_earth=dlon_earth-r360
if (dlon_earth<zero ) dlon_earth=dlon_earth+r360
dlat_earth = dlat_earth * deg2rad
dlon_earth = dlon_earth * deg2rad
if(regional)then
call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
if (outside) cycle
dlatmax=max(dlat,dlatmax)
dlonmax=max(dlon,dlonmax)
dlatmin=min(dlat,dlatmin)
dlonmin=min(dlon,dlonmin)
else
dlat = dlat_earth
dlon = dlon_earth
call grdcrd1(dlat,rlats,nlat,1)
call grdcrd1(dlon,rlons,nlon,1)
endif
! write(6,*)'dlat, dlon::',dlat,dlon
cdata_all(1,i)=0.0_r_kind ! time
cdata_all(2,i)=dlon ! grid relative longitude
cdata_all(3,i)=dlat ! grid relative latitude
cdata_all(4,i)=ptop_l(i) ! cloud top pressure (pa)
cdata_all(5,i)=teff_l(i) ! Cloud top temperature (K)
cdata_all(6,i)=phase_l(i) ! cloud cover
cdata_all(7,i)=lwp_l(i)
cdata_all(8,i)=usage
enddo
ilon=2
ilat=3
nchanl=0
nreal=maxdat
nread=numobs
ndata=numobs
nouse=0
call count_obs(numobs,maxdat,ilat,ilon,cdata_all,nobs)
write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
write(lunout) ((cdata_all(k,i),k=1,maxdat),i=1,numobs)
write(6,*)'NASA larcglb::',nreal,numobs
return
end subroutine read_NASA_LaRC_cloud
subroutine read_NASALaRC_cloud_bufr(satfile,atime,&
maxobs,numobs,ptop, teff, phase, lwp_iwp,lat, lon)
!
! PRGMMR: Ming Hu ORG: GSD DATE: 2010-07-09
!
! ABSTRACT:
! This routine read in NASA LaRC cloud products
! from a bufr file
!
! PROGRAM HISTORY LOG:
!
! variable list
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90 + EXTENSIONS
! MACHINE: wJET
!
!$$$
!
!_____________________________________________________________________
!
use kinds, only: r_kind,i_kind,r_double,r_single
implicit none
!
!
!
character(80):: hdstr='YEAR MNTH DAYS HOUR MINU SECO'
character(80):: obstr='CLATH CLONH CLDP HOCT CDTP EBBTH VILWC'
! CLDP | CLOUD PHASE
! HOCB | HEIGHT OF BASE OF CLOUD
! HOCT | HEIGHT OF TOP OF CLOUD (METERS)
! CDBP | PRESSURE AT BASE OF CLOUD
! CDTP | PRESSURE AT TOP OF CLOUD (PA)
! EBBTH | EQUIVALENT BLACK BODY TEMPERATURE (KELVIN)
! VILWC | VERTICALLY-INTEGRATED LIQUID WATER CONTENT
real(8) :: hdr(6),obs(7,1)
INTEGER(i_kind) :: ireadmg,ireadsb
character(8) subset
integer :: unit_in=10,idate,iret,nmsg,ntb
!
! For NASA LaRC
!
CHARACTER*40 satfile
INTEGER :: maxobs, numobs ! dimension
INTEGER(i_kind) :: obs_time
REAL*4 lat ( maxobs)
REAL*4 lon ( maxobs)
integer phase ( maxobs)
REAL*4 lwp_iwp ( maxobs)
REAL*4 teff ( maxobs)
REAL*4 ptop ( maxobs)
!
integer phase_tmp
!
!
! ** misc
! integer i,j,k
!
! integer :: status
character*10 atime
!
!**********************************************************************
!
open(24,file='NASA.bufrtable')
open(unit_in,file=trim(satfile),form='unformatted')
call openbf(unit_in,'IN',unit_in)
call dxdump(unit_in,24)
call datelen(10)
nmsg=0
ntb = 0
msg_report: do while (ireadmg(unit_in,subset,idate) == 0)
nmsg=nmsg+1
sb_report: do while (ireadsb(unit_in) == 0)
call ufbint(unit_in,hdr,6,1,iret,hdstr)
obs_time=int((hdr(1)-2000.0_r_kind)*100000000+hdr(2)*1000000+hdr(3)*10000+hdr(4)*100+hdr(5))
call ufbint(unit_in,obs,7,1,iret,obstr)
ntb = ntb+1
if(ntb > maxobs) then
write(*,*) 'Error: need to increase maxobs',maxobs, ntb
stop 1234
endif
lat(ntb)=obs(1,1)
lon(ntb)=obs(2,1)
phase_tmp=-9
if(obs(3,1) < 100.0_r_kind .and. obs(3,1) > -0.0_r_kind) phase_tmp=int(obs(3,1))
if (phase_tmp==0 .or. phase_tmp==1 .or. phase_tmp==2) then
if( (obs(5,1) < 1.e10_r_kind .and. obs(5,1) > 1.0_r_kind ) .and. &
(obs(6,1) < 1.e10_r_kind .and. obs(6,1) > 1.0_r_kind) ) then
phase(ntb) = 1 !cloud
lwp_iwp(ntb)=obs(7,1)
teff(ntb)=obs(6,1)
ptop(ntb)=obs(5,1)/100.0_r_kind ! pa to hPa
else
phase(ntb) = -9 ! bad data or no data
endif
elseif (phase_tmp==4) then
phase(ntb) = 0 !clear
ptop(ntb) = -20._r_single
else
phase(ntb) = -9 ! bad data or no data
endif
enddo sb_report
enddo msg_report
write(*,*) 'message/reports num=',nmsg,ntb
call closbf(unit_in)
numobs=ntb
write(atime,'(I10)') idate
end subroutine read_NASALaRC_cloud_bufr
subroutine read_NASALaRC_cloud_bufr_survey(satfile,east_time, west_time)
!
! PRGMMR: Ming Hu ORG: GSD DATE: 2010-07-09
!
! ABSTRACT:
! This routine read in NASA LaRC cloud products
! from a bufr file
!
! PROGRAM HISTORY LOG:
! 2015-05-12 s. liu - increase max_obstime from 10 to 20
!
! variable list
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90 + EXTENSIONS
! MACHINE: wJET
!
!$$$
!
!_____________________________________________________________________
!
use kinds, only: r_kind,i_kind,r_double
implicit none
!
!
!
character(80):: hdstr='YEAR MNTH DAYS HOUR MINU SECO'
real(8) :: hdr(6)
integer(i_kind) :: ireadmg,ireadsb
character(8) subset
integer(i_kind) :: unit_in=10,idate,iret,nmsg,ntb
!
! For NASA LaRC
!
CHARACTER*40, intent(in) :: satfile
integer(i_kind),intent(out) :: east_time, west_time
INTEGER(i_kind) :: obs_time
INTEGER(i_kind),parameter :: max_obstime=20
integer(i_kind) :: num_obstime_all(max_obstime)
integer(i_kind) :: num_subset_all(max_obstime)
integer(i_kind) :: num_obstime_hh(max_obstime)
integer(i_kind) :: num_obstime
!
integer :: i,ii,hhh
!
!**********************************************************************
!
num_obstime=0
hhh=99
open(24,file='NASA.bufrtable')
open(unit_in,file=trim(satfile),form='unformatted')
call openbf(unit_in,'IN',unit_in)
call dxdump(unit_in,24)
call datelen(10)
nmsg=0
msg_report: do while (ireadmg(unit_in,subset,idate) == 0)
ntb = 0
nmsg=nmsg+1
sb_report: do while (ireadsb(unit_in) == 0)
call ufbint(unit_in,hdr,6,1,iret,hdstr)
obs_time=int((hdr(1)-2000.0_r_kind)*100000000+hdr(2)*1000000+hdr(3)*10000+hdr(4)*100+hdr(5))
hhh=int(hdr(5))
ntb=ntb+1
enddo sb_report
! message inventory
if(num_obstime == 0 ) then
num_obstime=1
num_obstime_all(num_obstime)=obs_time
num_obstime_hh(num_obstime)=hhh
num_subset_all(num_obstime)= ntb
else
ii=0
DO i=1,num_obstime
if(num_obstime_all(i) == obs_time ) ii=i
ENDDO
if( ii > 0 .and. ii <=num_obstime) then
num_subset_all(ii)=num_subset_all(ii) + ntb
else
num_obstime=num_obstime+1
if(num_obstime> max_obstime) then
write(*,*) 'Error: too many message types'
write(*,*) 'Need to increase :max_obstime'
stop 1234
endif
num_obstime_all(num_obstime)=obs_time
num_obstime_hh(num_obstime)=hhh
num_subset_all(num_obstime)=num_subset_all(num_obstime)+ntb
endif
endif
enddo msg_report
write(*,*) 'message/reports num=',nmsg,ntb
call closbf(unit_in)
write(*,'(2x,a10,a10,a11)') 'time_level','subset_num'
DO i=1,num_obstime
write(*,'(i2,i12,i11,i10)') i,num_obstime_all(i),num_subset_all(i),num_obstime_hh(i)
ENDDO
! GOES EAST : 1815, 1845, 1915, 2045
! GOES WEST : 1830, 1900, 2030
east_time=0
west_time=0
DO i=1,num_obstime
if(num_subset_all(i) > 10) then
if(num_obstime_hh(i) == 15 .or. num_obstime_hh(i) == 45 ) then
if(east_time < num_obstime_all(i)) east_time=num_obstime_all(i)
endif
if(num_obstime_hh(i) == 30 .or. num_obstime_hh(i) == 0 ) then
if(west_time < num_obstime_all(i)) west_time=num_obstime_all(i)
endif
endif
ENDDO
write(*,*) 'east_time=',east_time
write(*,*) 'west_time=',west_time
end subroutine read_NASALaRC_cloud_bufr_survey