subroutine read_goesimgr_skycover(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis,& prsl_full,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_goesimgr_skycover read GOES Imager sky cover product ! prgmmr: Jacob Carley date: 2014-11-07 ! ! abstract: This routine reads GOES Imager sky cover data from bufr_d files. ! It also has options to thin the data by using conventional ! thinning programs, though data are only thinned in 2D and not 3D ! since sky cover is a 2D data set. ! When running the gsi in regional mode, the code only ! retains those observations that fall within the regional ! domain ! ! program history log: ! 2014-11-07 J. Carley - Initial code ! 2015-03-06 C. Thomas - Add thin4d logical for removal of time thinning ! 2016-01-11 D. Keyser - Enable use of efclam dump as the primary observation ! source file. However, if ob is missing, will look for it in old ! BUFR mnemonic sequence ! 2016-03-11 j. guo - Fixed {dlat,dlon}_earth_deg in the obs data stream ! 2016-04-22 M. Pondeca - Replace "if (goescld(3)==bmiss)" condition with "if (goescld(3) > r0_01_bmiss)" ! 2019-06-17 M. Morris - Update adjust_goescldobs to reject clear cloud obs over water at night (tcamt_qc==8) ! 2019-12-05 M. Morris - Update adjust_goescldobs to reject ALL clear cloud cover obs at night (tcamt_qc==8) ! ! input argument list: ! ithin - flag to thin data ! rmesh - thinning mesh size (km) ! gstime - analysis time in minutes from reference date ! infile - unit from which to read BUFR data ! lunout - unit to which to write data for further processing ! obstype - observation type to process ! twind - input group time window (hours) ! sis - satellite/instrument/sensor indicator ! ! output argument list: ! nread - number of obs read ! ndata - number of obs retained for further processing ! nodata - number of obs retained for further processing ! nobs - array of observations on each subdomain for each processor ! ! attributes: ! language: f95/2003 ! machine: WCOSS ! !$$$ use kinds, only: r_single,r_kind,r_double,i_kind use constants, only: zero,one_tenth,one,deg2rad,half,& three,four, r60inv,r10,r100,r2000 use convinfo, only: nconvtype, & icuse,ictype,ioctype,& ithin_conv,rmesh_conv,pmesh_conv,ctwind use convthin, only: make3grids,map3grids,del3grids,use_all use gridmod, only: regional,nlon,nlat,nsig,tll2xy,txy2ll,& rlats,rlons use deter_sfc_mod, only: deter_sfc2 use obsmod, only: bmiss,ran01dom use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen,time_4dvar,thin4d use adjust_cloudobs_mod, only: adjust_goescldobs use mpimod, only: npe implicit none ! Declare passed variables character(len=*) ,intent(in ) :: infile,obstype character(len=20) ,intent(in ) :: sis integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: nread,ndata,nodata integer(i_kind),dimension(npe) ,intent(inout) :: nobs real(r_kind) ,intent(in ) :: twind,gstime real(r_kind),dimension(nlat,nlon,nsig),intent(in ) :: prsl_full ! Declare local parameters real(r_kind),parameter:: r90 = 90.0_r_kind real(r_kind),parameter:: r0_01 = 0.01_r_kind real(r_kind),parameter:: r0_1_bmiss=one_tenth*bmiss real(r_kind),parameter:: r0_01_bmiss=r0_01*bmiss real(r_kind),parameter:: r1200= 1200.0_r_kind real(r_kind),parameter:: r6= 6.0_r_kind real(r_kind),parameter:: r360 = 360.0_r_kind character(8),parameter:: cspval= '88888888' ! Declare local variables character(len=80) :: hdrstr,goescldstr,goescldstr_new character(len=8) :: subset character(len=22) :: myname character(len=8) :: c_prvstg,c_sprvstg ,c_station_id integer(i_kind) :: nmsub,ireadmg,ireadsb,nreal,nc,i,lunin,nmsg,ntb integer(i_kind) :: iret,kx,pflag,nlevp,nmind,levs,idomsfc integer(i_kind) :: low_cldamt_qc,mid_cldamt_qc,hig_cldamt_qc,tcamt_qc integer(i_kind) :: ithin,klat1,klon1,klonp1,klatp1,kk,k,ilat,ilon,nchanl integer(i_kind) :: iout,ntmp,iiout,maxobs,icount,itx,iuse,idate,ierr integer(i_kind),dimension(5) :: idate5 integer(i_kind),allocatable,dimension(:):: isort,iloc real(r_kind) :: dlat,dlon,dlat_earth,dlon_earth,toff,t4dv real(r_kind) :: dlat_earth_deg,dlon_earth_deg real(r_kind) :: dx,dx1,dy,dy1,w00,w10,w01,w11,crit1,timedif,tdiff real(r_kind) :: rmesh,pmesh,xmesh,tcamt,tcamt_oe,ff10,tsavg real(r_kind) :: rminobs,ppb real(r_kind) :: low_cldamt,mid_cldamt,hig_cldamt,usage,zz,sfcr,rstation_id real(r_kind),allocatable,dimension(:):: presl_thin real(r_kind),dimension(nsig):: presl real(r_kind),allocatable,dimension(:,:):: cdata_all,cdata_out real(r_double),dimension(9):: hdr real(r_double),dimension(3):: goescld logical :: outside,ithinp,luse real(r_double),dimension(1,1):: r_prvstg,r_sprvstg ! equivalence to handle character names equivalence(rstation_id,c_station_id) lunin=11_i_kind ithin=-9_i_kind rmesh=-99.999_r_kind myname='READ_GOESIMGR_SKYCOVER' hdrstr='NUL YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH' goescldstr='SAID TOCC TOCC_AVG' goescldstr_new='SAID ECAS ECAM' nreal=20 nc=0 conv: do i=1,nconvtype if(trim(obstype) == trim(ioctype(i)) .and. ictype(i)==154_i_kind) then nc=i exit conv end if end do conv if(nc == 0)then write(6,*) myname,' no matching obstype found in convinfo ',obstype return end if ! Try opening bufr file, if unable print error to the screen ! and return to read_obs.F90 open(lunin,file=trim(infile),form='unformatted',iostat=ierr) if (ierr/=0) then write(6,'(A)')myname,':ERROR: Trouble opening input file: ',trim(infile),' returning to read_obs.F90...' return end if call openbf(lunin,'IN',lunin) call datelen(10) ! Set up thinning params use_all = .true. ithin=ithin_conv(nc) if (ithin > 0 ) then rmesh=rmesh_conv(nc) pmesh=pmesh_conv(nc) use_all = .false. if(pmesh > zero) then pflag=1 nlevp=r1200/pmesh else pflag=0 nlevp=nsig endif pflag=0 nlevp=nsig xmesh=rmesh call make3grids(xmesh,nlevp) if (.not.use_all) then allocate(presl_thin(nlevp)) if (pflag==1) then do k=1,nlevp presl_thin(k)=(r1200-(k-1)*pmesh)*one_tenth enddo endif endif write(6,'(A,1x,A,1x,A,I4,1x,f8.2,1x,I2,1x,I3,1x,f8.2,1x,I3)')myname,': ioctype(nc),ictype(nc),rmesh,pflag,nlevp,pmesh,nc ',& trim(ioctype(nc)),ictype(nc),rmesh,pflag,nlevp,pmesh,nc endif nmsg=0 ntb = 0 ! Find number of reports and messages so we know how much memory to allocate do while (ireadmg(lunin,subset,idate) == 0) ! Time offset if(nmsg == 0) call time_4dvar(idate,toff) nmsg=nmsg+1 ntb = ntb + nmsub(lunin) !nmsub is a bufrlib function which returns the number of subsets in ! a bufr message open for input via a previous call to a bufrlib ! routine readmg or equivalent. The subsets are not required to be read (saves time). end do maxobs=ntb allocate(cdata_all(nreal,maxobs),isort(maxobs)) isort = 0 cdata_all=zero nread=0 nchanl=0 ilon=2 ilat=3 ntb=0 call closbf(lunin) close(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) loop_msg: do while (ireadmg(lunin,subset,idate) == 0) loop_readsb: do while (ireadsb(lunin) == 0) ntb=ntb+1 ! - Extract type, date, and location information call ufbint(lunin,hdr,9,1,iret,hdrstr) ! - Compare relative obs time with window. If obs ! - falls outside of window, don't use this obs idate5(1) = hdr(2) ! year idate5(2) = hdr(3) ! month idate5(3) = hdr(4) ! day idate5(4) = hdr(5) ! hours idate5(5) = hdr(6) ! minutes call w3fs21(idate5,nmind) rminobs=real(nmind,8)+(real(hdr(7),8)*r60inv)!convert the seconds of the ob to minutes and store to rminobs t4dv = (rminobs-real(iwinbgn,r_kind))*r60inv tdiff=(rminobs-gstime)*r60inv !GS time is the analysis time in minutes from w3fs21 if (l4dvar.or.l4densvar) then if (t4dvwinlen) cycle loop_readsb else ! - Check to make sure ob is within convinfo time window (ctwind) and ! - is within overwall time window twind (usually +-3) if( (abs(tdiff) > ctwind(nc)) .or. (abs(tdiff) > twind) )cycle loop_readsb endif kx=999_i_kind !hardwire typ to 999 if(abs(hdr(8))>r90 .or. abs(hdr(9))>r360) cycle loop_readsb if(hdr(9)== r360)hdr(9)=hdr(9)-r360 if(hdr(9) < zero)hdr(9)=hdr(9)+r360 dlon_earth_deg = hdr(9) dlat_earth_deg = hdr(8) dlon_earth=hdr(9)*deg2rad dlat_earth=hdr(8)*deg2rad nread=nread+1 if(regional)then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) ! convert to rotated coordinate if(outside) cycle loop_readsb ! check to see if outside regional domain else dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif ! Read in the obs goescld=bmiss call ufbint(lunin,goescld,3,1,levs,goescldstr_new) if (goescld(3) > r0_01_bmiss) then ! if ob is missing, look for it in old BUFR mnemonic sequence goescld=bmiss call ufbint(lunin,goescld,3,1,levs,goescldstr) if (goescld(3) > r0_01_bmiss) cycle loop_readsb !If obs are missing, cycle endif c_prvstg=cspval c_sprvstg=cspval ! - Set station ID rstation_id=goescld(1) ithin=ithin_conv(nc) ithinp = ithin > 0 .and. pflag /= 0 ! - Thin in vertical - note we can only thin in the horizontal ! - since sky cover is a 2D field. So this branch should never run ! - unless we get info about the vertical location of the clouds in the ! - future. Leaving here as a 'just-in-case' measure. if(ithinp )then ! Interpolate guess pressure profile to observation location klon1= int(dlon); klat1= int(dlat) dx = dlon-klon1; dy = dlat-klat1 dx1 = one-dx; dy1 = one-dy w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon) if (klon1==0) klon1=nlon klatp1=min(nlat,klat1+1); klonp1=klon1+1 if (klonp1==nlon+1) klonp1=1 do kk=1,nsig presl(kk)=w00*prsl_full(klat1 ,klon1 ,kk) + & w10*prsl_full(klatp1,klon1 ,kk) + & w01*prsl_full(klat1 ,klonp1,kk) + & w11*prsl_full(klatp1,klonp1,kk) end do end if iuse=icuse(nc) ! General block for data thinning - if requested if (ithin > 0 .and. iuse >=0) then ntmp=ndata ! counting moved to map3gridS ! - Set data quality index for thinning if (thin4d) then timedif = zero else timedif=abs(t4dv-toff) endif crit1 = timedif/r6+half ! - simple 1-to-1 mapping of vertical levels when no thinning in the vertical if (pflag==0) then do kk=1,nsig presl_thin(kk)=presl(kk) end do endif ppb=one_tenth*1013.25_r_kind !number is irrelevant for 2D - set to standard SLP -> 1013.25 and convert from mb to cb call map3grids(-1,pflag,presl_thin,nlevp,dlat_earth,dlon_earth,& ppb,crit1,ndata,iout,ntb,iiout,luse,.false.,.false.) if (.not. luse) cycle loop_readsb if(iiout > 0) isort(iiout)=0 if (ndata > ntmp) then nodata=nodata+1 endif isort(ntb)=iout else ! - no thinnning ndata=ndata+1 nodata=nodata+1 iout=ndata isort(ntb)=iout endif !- Set usage variable usage = 0 if(iuse <= 0)usage=r100 ! Get information from surface file necessary for conventional data here call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz) ! - Obtain the ob and tune the QC marks for ob error tuning a bit later call adjust_goescldobs(goescld(3),tdiff,dlat_earth,dlon_earth, & low_cldamt,low_cldamt_qc,mid_cldamt,mid_cldamt_qc, & hig_cldamt,hig_cldamt_qc,tcamt,tcamt_qc) if(tcamt_qc==15 .or. tcamt_qc==12 .or. tcamt_qc==9 .or. tcamt_qc==8) usage=r100 tcamt_oe=20.0_r_kind if(tcamt_qc==1) tcamt_oe=tcamt_oe*1.25_r_kind if(tcamt_qc==2) tcamt_oe=tcamt_oe*1.50_r_kind if(tcamt_qc==3) tcamt_oe=tcamt_oe*1.75_r_kind cdata_all( 1,iout)=tcamt_oe ! obs error cdata_all( 2,iout)=dlon ! grid relative longitude cdata_all( 3,iout)=dlat ! grid relative latitude cdata_all( 4,iout)=tcamt ! total cloud amount (%) cdata_all( 5,iout)=rstation_id ! station ID cdata_all( 6,iout)=t4dv ! time cdata_all( 7,iout)=nc ! type cdata_all( 8,iout)=tcamt_qc ! quality mark cdata_all( 9,iout)=usage ! usage parameter cdata_all(10,iout)=idomsfc ! dominate surface type cdata_all(11,iout)=tsavg ! skin temperature cdata_all(12,iout)=ff10 ! 10 meter wind factor cdata_all(13,iout)=sfcr ! surface roughness cdata_all(14,iout)=dlon_earth_deg ! earth relative longitude (degrees) cdata_all(15,iout)=dlat_earth_deg ! earth relative latitude (degrees) cdata_all(16,iout)=bmiss ! station elevation (m) cdata_all(17,iout)=bmiss ! observation height (m) cdata_all(18,iout)=zz ! terrain height at ob location cdata_all(19,iout)=r_prvstg(1,1) ! provider name cdata_all(20,iout)=r_sprvstg(1,1) ! subprovider name enddo loop_readsb enddo loop_msg ! Close unit to bufr file call closbf(lunin) ! Deallocate arrays used for thinning data if (.not.use_all) then deallocate(presl_thin) call del3grids endif ! Normal exit ! Write header record and data to output file for further processing allocate(iloc(ndata)) icount=0 do i=1,maxobs if(isort(i) > 0)then icount=icount+1 iloc(icount)=isort(i) end if end do if(ndata /= icount)then write(6,*) myname,': ndata and icount do not match STOPPING...ndata,icount ',ndata,icount call stop2(50) end if allocate(cdata_out(nreal,ndata)) do i=1,ndata itx=iloc(i) do k=1,nreal cdata_out(k,i)=cdata_all(k,itx) end do end do deallocate(iloc,isort,cdata_all) call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata write(lunout) cdata_out deallocate(cdata_out) if (ndata == 0) then write(6,*)myname,'no read_goesimgr_skycover data' endif close(lunin) end subroutine read_goesimgr_skycover