program process_NASALaRC_cloud
!
! PRGMMR: Ming Hu ORG: GSD DATE: 2009-09-04
!
! ABSTRACT:
! This routine read in NASA LaRC cloud products and
! interpolate them into GSI mass grid
!
! PROGRAM HISTORY LOG:
!
! variable list
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90 + EXTENSIONS
! MACHINE: wJET
!
!$$$
!
!_____________________________________________________________________
!
use mpi
use kinds, only: r_kind,i_kind,r_single
use map_utils
use misc_definitions_module , only : PROJ_LC, PROJ_ROTLL
use constants_module ,only : EARTH_RADIUS_M
use constants, only: init_constants_derived, deg2rad
use gridmod_gsimap ,only : nlon,nlat,init_general_transform,tll2xy,txy2ll
implicit none
!
INCLUDE 'netcdf.inc'
!
! MPI variables
integer :: npe, mype,ierror
!
real :: rad2deg = 180.0/3.1415926
!
character*256 output_file
!
! grid
! integer(i_kind) :: nlon,nlat
real,allocatable:: xlon(:,:) !
real,allocatable:: ylat(:,:) !
real(r_kind),allocatable:: rxlon(:,:) !
real(r_kind),allocatable:: rylat(:,:) !
real :: userDX, userDY, CEN_LAT, CEN_LON
real :: userTRUELAT1,userTRUELAT2,MOAD_CEN_LAT,STAND_LON
integer :: MAP_PROJ
type (proj_info) :: proj_stack
REAL :: truelat1, truelat2, stdlon, lat1, lon1, r_earth
REAL :: knowni, knownj, dx
REAL :: user_known_x,user_known_y
CHARACTER*180 geofile
!
! For NASA LaRC
!
CHARACTER*180 workPath
CHARACTER*80 satfile
INTEGER :: nxp, nyp ! dimension
! ****VARIABLES FOR THIS NETCDF FILE****
!
CHARACTER*24 :: cbase_time
INTEGER(i_kind) :: base_time
INTEGER(i_kind) :: ibase_year,ibase_month,ibase_day,ibase_hour,ihour
INTEGER(i_kind) :: icycle_year,icycle_month,icycle_day,icycle_hour
REAL*8 time_offset
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
!
REAL(r_single), allocatable :: w_pcld(:,:)
REAL(r_single), allocatable :: w_tcld(:,:)
REAL(r_single), allocatable :: w_frac(:,:)
REAL(r_single), allocatable :: w_lwp (:,:)
integer(i_kind),allocatable :: nlev_cld(:,:)
!
!
! Working
integer nfov
parameter (nfov=160)
real, allocatable :: Pxx(:,:,:),Txx(:,:,:),WPxx(:,:,:)
real,allocatable :: xdist(:,:,:), xxxdist(:)
real fr,sqrt, qc, type
integer,allocatable :: PHxx(:,:,:),index(:,:), jndex(:)
integer ixx,ii,jj,med_pt,igrid,jgrid &
,ncount,ncount1,ncount2,ii1,jj1,nobs,n
!
! namelist
!
integer :: analysis_time
integer :: ioption
character(len=100) :: bufrfile
integer(i_kind) :: npts_rad
namelist/setup/ analysis_time, ioption, npts_rad,bufrfile
!
!
! ** misc
real(r_kind) :: xc ! x-grid coordinate (grid units)
real(r_kind) :: yc ! y-grid coordinate (grid units)
real(r_kind) :: rlon ! earth longitude (radians)
real(r_kind) :: rlat ! earth latitude (radians)
logical ::outside ! .false., then point is inside x-y domain
! .true., then point is outside x-y domain
integer i,j,k,ipt,jpt,cfov
Integer nf_status,nf_fid,nf_vid
integer :: NCID
integer(i_kind) :: east_time, west_time
integer :: isat, maxobs,numobs
integer :: status
character*10 atime
logical :: ifexist
!**********************************************************************
!
! END OF DECLARATIONS....start of program
!
! MPI setup
call MPI_INIT(ierror)
call MPI_COMM_SIZE(mpi_comm_world,npe,ierror)
call MPI_COMM_RANK(mpi_comm_world,mype,ierror)
if(mype==0) then
! get namelist
!
analysis_time=2018051718
bufrfile='NASALaRCCloudInGSI.bufr'
npts_rad=1
! * ioption = 1 is nearest neighrhood
! * ioption = 2 is median of cloudy fov
ioption = 2
inquire(file='namelist_nasalarc', EXIST=ifexist )
if(ifexist) then
open(10,file='namelist_nasalarc',status='old')
read(10,setup)
close(10)
write(*,*) 'Namelist setup are:'
write(*,setup)
else
write(*,*) 'No namelist file exist, use default values'
write(*,*) "analysis_time,bufrfile,npts_rad,ioption"
write(*,*) analysis_time, trim(bufrfile),npts_rad,ioption
endif
! set geogrid fle name
!
call init_constants_derived
workPath='./'
write(geofile,'(a,a)') trim(workPath), 'geo_em.d01.nc'
write(*,*) 'geofile', trim(geofile)
call GET_DIM_ATT_geo(geofile,NLON,NLAT)
write(*,*) 'NLON,NLAT',NLON,NLAT
call GET_MAP_ATT_geo(geofile, userDX, userDY, CEN_LAT, CEN_LON, &
userTRUELAT1,userTRUELAT2,MOAD_CEN_LAT,STAND_LON,MAP_PROJ)
write(*,*) userDX, userDY, CEN_LAT, CEN_LON
write(*,*) userTRUELAT1,userTRUELAT2,MOAD_CEN_LAT,STAND_LON,MAP_PROJ
!
! get GSI horizontal grid in latitude and longitude
!
allocate(xlon(nlon,nlat),rxlon(nlon,nlat))
allocate(ylat(nlon,nlat),rylat(nlon,nlat))
call OPEN_geo(geofile, NCID)
call GET_geo_sngl_geo(NCID,Nlon,Nlat,ylat,xlon)
call CLOSE_geo(NCID)
!
! setup map
!
! if (MAP_PROJ == PROJ_LC) then
! user_known_x = (NLON+1)/2.0
! user_known_y = (NLAT+1)/2.0
! call map_init(proj_stack)
!
! call map_set(MAP_PROJ, proj_stack, &
! truelat1=userTRUELAT1, &
! truelat2=userTRUELAT2, &
! stdlon=STAND_LON, &
! lat1=CEN_LAT, &
! lon1=CEN_LON, &
! knowni=user_known_x, &
! knownj=user_known_y, &
! dx=userDX, &
! r_earth=earth_radius_m)
! else
mype=0
rylat=ylat*deg2rad
rxlon=xlon*deg2rad
call init_general_transform(rylat,rxlon,mype)
! endif
!
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
! maxobs=(1800*700 + 1500*850)*1
satfile='NASA_LaRC_cloud.bufr'
call read_NASALaRC_cloud_bufr_survey(satfile,east_time, west_time,maxobs)
allocate(lat_l(maxobs))
allocate(lon_l(maxobs))
allocate(ptop_l(maxobs))
allocate(teff_l(maxobs))
allocate(phase_l(maxobs))
allocate(lwp_l(maxobs))
ptop_l=-9.
teff_l=-9.
lat_l =-9.
lon_l =-9.
lwp_l =-9.
phase_l=-9
call read_NASALaRC_cloud_bufr(satfile,atime,east_time, west_time, &
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
if (phase_l(i).eq.4) phase_l(i) = 0 ! clear
if (phase_l(i).eq.5) phase_l(i) = -9 ! bad data
! equivalent to "no data"
! if (phase_l(i).eq.-9) ptop_l(i) = -9.
if (phase_l(i).eq.0) ptop_l(i) = -20.
enddo
! -----------------------------------------------------------
! -----------------------------------------------------------
! Map each FOV onto RR grid points
! -----------------------------------------------------------
! -----------------------------------------------------------
write(*,*) 'The number of impact point is=',npts_rad
do ipt=1,numobs
if (phase_l(ipt).ge.0) then
! Indicates there is some data (not missing)
! if (MAP_PROJ == PROJ_LC) then
! call latlon_to_ij(proj_stack, 90.0-lat_l(ipt,jpt), lon_l(ipt,jpt), xc, yc)
! else
rlon=lon_l(ipt)*deg2rad
!mhu rlat=(90.0-lat_l(ipt,jpt))*deg2rad
rlat=lat_l(ipt)*deg2rad
call tll2xy(rlon,rlat,xc,yc)
! call txy2ll(xc,yc,rlon,rlat)
! endif
! * Compute RR grid x/y at lat/lon of cloud data
! -----------------------------------------------------------
! * XC,YC should be within RR boundary, i.e., XC,YC >0
ii1 = int(xc+0.5)
jj1 = int(yc+0.5)
do jj = max(1,jj1-npts_rad), min(nlat,jj1+npts_rad)
if (jj1-1.ge.1 .and. jj1+1.le.nlat) then
do ii = max(1,ii1-npts_rad), min(nlon,ii1+npts_rad)
if (ii1-1.ge.1 .and. ii1+1.le.nlon) then
! if(XC .ge. 1. .and. XC .lt. nlon .and. &
! YC .ge. 1. .and. YC .lt. nlat) then
! ii1 = int(xc+0.5)
! jj1 = int(yc+0.5)
! ii=ii1
! jj=jj1
! * We check multiple data within gridbox
if (index(ii,jj).lt.nfov) then
index(ii,jj) = index(ii,jj) + 1
Pxx(ii,jj,index(ii,jj)) = Ptop_l(ipt)
Txx(ii,jj,index(ii,jj)) = Teff_l(ipt)
PHxx(ii,jj,index(ii,jj)) = phase_l(ipt)
WPxx(ii,jj,index(ii,jj)) = lwp_l(ipt)
xdist(ii,jj,index(ii,jj)) = sqrt( &
(XC+1-ii)**2 + (YC+1-jj)**2)
else
write(6,*) 'ALERT: too many data in one grid, increase nfov'
write(6,*) nfov, ii,jj
end if
endif
enddo ! ii
endif
enddo ! jj
! endif ! observation is in the domain
endif ! phase_l >= 0
enddo ! ipt
deallocate(lat_l,lon_l,ptop_l,teff_l,phase_l,lwp_l)
!
write(6,*) 'The max index number is: ', maxval(index)
allocate(w_pcld(nlon,nlat))
allocate(w_tcld(nlon,nlat))
allocate(w_frac(nlon,nlat))
allocate(w_lwp(nlon,nlat))
allocate(nlev_cld(nlon,nlat))
w_pcld=99999.
w_tcld=99999.
w_frac=99999.
w_lwp=99999.
nlev_cld = 99999
do jj = 1,nlat
do ii = 1,nlon
if ((index(ii,jj) .ge. 1 .and. index(ii,jj) .lt. 3) .and. userDX < 7000.0 ) then
w_pcld(ii,jj) = Pxx(ii,jj,1) ! hPa
w_tcld(ii,jj) = Txx(ii,jj,1) ! K
w_lwp(ii,jj) = WPxx(ii,jj,1) ! g/m^2
w_frac(ii,jj) = 1
nlev_cld(ii,jj) = 1
if (w_pcld(ii,jj).eq.-20) then
w_pcld(ii,jj) = 1013. ! hPa - no cloud
w_frac(ii,jj)=0.0
nlev_cld(ii,jj) = 0
end if
elseif(index(ii,jj) .ge. 3) then
! * We decided to use nearest neighborhood for ECA values,
! * a kind of convective signal from GOES platform...
!
! * 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 sortmed(xxxdist,index(ii,jj),jndex,fr)
w_pcld(ii,jj) = Pxx(ii,jj,jndex(1))
w_tcld(ii,jj) = Txx(ii,jj,jndex(1))
w_lwp(ii,jj) = WPxx(ii,jj,jndex(1))
endif
! * Sort to find median value
if(ioption .eq. 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
w_pcld(ii,jj) = Pxx(ii,jj,jndex(med_pt)) ! hPa
w_tcld(ii,jj) = Txx(ii,jj,jndex(med_pt)) ! K
w_lwp(ii,jj) = WPxx(ii,jj,jndex(med_pt)) ! g/m^2
endif ! pick median
! missing pcld
if (w_pcld(ii,jj).eq.-20) then
w_pcld(ii,jj) = 1013. ! hPa - no cloud
w_frac(ii,jj)=0.0
nlev_cld(ii,jj) = 0
! cloud fraction based on phase (0 are clear), what about -9 ????
elseif( w_pcld(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
endif ! index > 3
enddo !ii
enddo !jj
write (6,*) 'index'
write (6,'(10i10)') (index(300,jj),jj=1,nlat,10)
write (6,*) 'w_pcld'
write (6,'(10f10.2)') (w_pcld(300,jj),jj=1,nlat,10)
write (6,*) 'w_tcld'
write (6,'(10f10.2)') (w_tcld(300,jj),jj=1,nlat,10)
write (6,*) 'w_frac'
write (6,'(10f10.2)') (w_frac(300,jj),jj=1,nlat,10)
write (6,*) 'w_lwp '
write (6,'(10f10.2)') (w_lwp (300,jj),jj=1,nlat,10)
write (6,*) 'nlev_cld'
write (6,'(10I10)') (nlev_cld(300,jj),jj=1,nlat,10)
OPEN(15, file='NASALaRC_cloud.bin',form='unformatted')
write(15) nlon,nlat
write(15) xlon
write(15) ylat
write(15) index
write(15) w_pcld
write(15) w_tcld
write(15) w_frac
write(15) w_lwp
write(15) nlev_cld
CLOSE(15)
!
! write out results
!
call write_bufr_NASALaRC(bufrfile,analysis_time,nlon,nlat,userDX,index,w_pcld,w_tcld,w_frac,w_lwp,nlev_cld)
!
write(6,*) "=== RAPHRRR PREPROCCESS SUCCESS ==="
endif ! if mype==0
call MPI_FINALIZE(ierror)
! DO j=1,nlat
! DO I=1,nlon
! if( w_lwp(i,j) > 99998.0 ) w_lwp(i,j) = 0.0
! ENDDO
! ENDDO
!
end program process_NASALaRC_cloud
subroutine read_NASALaRC_cloud_bufr(satfile,atime,east_time, west_time, &
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
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 :: ireadmg,ireadsb
character(8) subset
integer :: unit_in=10,idate,iret,nmsg,ntb
!
! For NASA LaRC
!
CHARACTER*40 satfile
integer(i_kind) :: east_time, west_time
INTEGER :: maxobs, numobs ! dimension
INTEGER(i_kind) :: obs_time
REAL*8 time_offset
REAL*4 lat ( maxobs)
REAL*4 lon ( maxobs)
integer phase ( maxobs)
REAL*4 lwp_iwp ( maxobs)
REAL*4 teff ( maxobs)
REAL*4 ptop ( maxobs)
!
!
! ** misc
! integer i,j,k
! Integer nf_status,nf_fid,nx,ny,nf_vid
!
! 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)*100000000+hdr(2)*1000000+hdr(3)*10000+hdr(4)*100+hdr(5))
call ufbint(unit_in,obs,7,1,iret,obstr)
if(obs_time == east_time .or. obs_time == west_time ) then
if(abs(obs(3,1) -4.0) < 1.e-4) then
obs(7,1)=99999. ! clear
obs(6,1)=99999. ! clear
obs(5,1)=101300.0 ! clear (hpa)
endif
if(obs(5,1) < 1.e7 .and. obs(5,1) > 100.0 ) then
if(obs(6,1) < 1.e7 .and. obs(6,1) > 10.0) then
ntb = ntb+1
if(ntb > maxobs) then
write(*,*) 'ALERT: need to increase maxobs',maxobs, ntb
ntb = maxobs
endif
lwp_iwp(ntb)=99999.0
lat(ntb)=99999.0
lon(ntb)=99999.0
phase(ntb)=99999
teff(ntb)=99999.0
ptop(ntb)=99999.0
if(obs(1,1) < 1.e9) lat(ntb)=real(obs(1,1))
if(obs(2,1) < 1.e9) lon(ntb)=real(obs(2,1))
if(obs(3,1) < 1.e9) phase(ntb)=int(obs(3,1))
if(obs(7,1) < 1.e9) lwp_iwp(ntb)=real(obs(7,1))
if(obs(6,1) < 1.e9) teff(ntb)=real(obs(6,1))
if(obs(5,1) < 1.e9) ptop(ntb)=real(obs(5,1))/100.0 ! pa to hpa
endif
endif
endif ! east_time, west_time
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 sortmed(p,n,is)
real p(n)
integer is(n)
! * count cloudy fov
real f
integer cfov
cfov = 0
do i=1,n
! - changed for NASA LaRC, p set = -9 for clear FOVs
if(p(i) .gt. 0.) cfov = cfov + 1
enddo
f = float(cfov)/(max(1,n))
! cloud-top pressure is sorted high cld to clear
nm1 = n-1
do 10 i=1,nm1
ip1 = i+1
do 10 j=ip1,n
if(p(i).le.p(j)) goto 10
temp = p(i)
p(i) = p(j)
p(j) = temp
iold = is(i)
is(i) = is(j)
is(j) = iold
10 continue
return
end subroutine sortmed
subroutine read_NASALaRC_cloud_bufr_survey(satfile,east_time, west_time,maxobs)
!
! 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
implicit none
!
!
!
character(80):: hdstr='YEAR MNTH DAYS HOUR MINU SECO'
real(8) :: hdr(6)
INTEGER :: ireadmg,ireadsb
character(8) subset
integer :: 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,intent(out) :: maxobs
INTEGER :: numobs ! dimension
INTEGER(i_kind) :: obs_time
REAL*8 time_offset
INTEGER(i_kind),parameter :: max_obstime=30
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
!
character*10 atime
integer :: i,ii,hhh
integer :: numobs_east, numobs_west
!
!**********************************************************************
!
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)*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)) then
east_time=num_obstime_all(i)
numobs_east=num_subset_all(i)
endif
endif
if(num_obstime_hh(i) == 30 .or. num_obstime_hh(i) == 0 ) then
if(west_time < num_obstime_all(i)) then
west_time=num_obstime_all(i)
numobs_west=num_subset_all(i)
endif
endif
endif
ENDDO
write(*,*) 'east_time and number=',east_time,numobs_east
write(*,*) 'west_time and number=',west_time,numobs_west
maxobs=numobs_west+numobs_east
maxobs=maxobs+int(maxobs*0.2)
write(*,*) 'maxobs=',maxobs
end subroutine read_NASALaRC_cloud_bufr_survey