subroutine read_goesglm(nread,ndata,nodata,infile,obstype,lunout,twindin,sis) !$$$ subprogram documentation block ! . . . . ! subprogram: read_goesglm reads lightning obs from a BUFR file ! prgmmr: k apodaca ! org: CSU/CIRA, Data Assimilation Group ! date: 2015-03-12 ! ! abstract: This routine reads lightning data (Earth-relative location, frequency) from ! file and prepares it for assimiliation in the form of lightning flash rate ! (#hits km-2 hr-1). ! ! Note: when running GSI in regional mode, the code only ! retains those observations that fall within the regional ! domain ! ! program history log: ! 2015-06-12 zupanski - include the convert_to_flash_rate subroutine to convert ! lightning strike observations into lightning flash rate ! (#hits km-2 hr-1). ! 2015-11-18 apodaca - include the convert_time subroutine to deal with computer ! precision dependencies. ! 2018-02-07 apodaca - add further documentation ! ! input argument list: ! infile - unit from which to read BUFR data ! obstype - observation type to process ! lunout - unit to which to write data for further processing ! ! output argument list: ! nread - number of type "obstype" observations read ! nodata - number of individual "obstype" observations read ! ndata - number of type "obstype" observations retained for further processing ! twindin - input group time window (hours) ! sis - satellite/instrument/sensor indicator ! ! attributes: ! language: Fortran 90 and/or above ! machine: ! !$$$ use kinds, only: r_single,r_kind,r_double,i_kind use constants, only: zero,one_tenth,one,deg2rad,& three,rad2deg,& r60inv,ten use gridmod, only: diagnostic_reg,wrf_mass_regional,regional,nlon,nlat,& tll2xy,txy2ll,& rlats,rlons use lightinfo, only: iuse_light,nlighttype use obsmod, only: iadate use obsmod, only: offtime_data use gsi_4dvar, only: l4dvar,l4densvar,time_4dvar,winlen implicit none ! Declare passed variables character(len=*) ,intent(in ) :: infile,obstype character(len=*) ,intent(in ) :: sis integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: nread,ndata real(r_kind) ,intent(in ) :: twindin integer(i_kind) :: nodata ! Declare local parameters real(r_kind),parameter:: r6 = 6.0_r_kind real(r_kind),parameter:: r90 = 90.0_r_kind real(r_kind),parameter:: r180 = 180.0_r_kind real(r_kind),parameter:: r360 = 360.0_r_kind !--- Declare local variables logical lob logical outside character(40) hdstr,oestr,qcstr character(10) date character(8) subset character(1) sidchr(8) integer(i_kind) ireadmg,ireadsb,icntpnt,icount integer(i_kind) lunin,i integer(i_kind) itx integer(i_kind) ihh,idd,idate,iret,im,iy,k integer(i_kind) nchanl,nreal,ilat,ilon integer(i_kind) lqm integer(i_kind) iout integer(i_kind) ntest,nvtest integer(i_kind) minobs,minan integer(i_kind) ntb integer(i_kind) nmsg ! message index integer(i_kind),parameter :: maxobs=2000000 integer(i_kind),dimension(5):: idate5 integer(i_kind),allocatable,dimension(:):: isort,iloc real(r_kind) time real(r_kind) usage real(r_kind) loe,lmerr real(r_kind) time_correction real(r_kind) dlat,dlon,dlat_earth,dlon_earth real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00 real(r_kind) vdisterrmax real(r_kind) timex,timeobs,toff,t4dv,zeps real(r_kind),allocatable,dimension(:,:):: cdata_all,cdata_out !--- flash rate real(r_kind),allocatable,dimension(:,:):: cdata_flash,cdata_flash_h integer(i_kind) :: ndata_flash,ndata_flash_h real(r_double) rstation_id real(r_double),dimension(3):: hdr real(r_double),dimension(1,1):: qcmark,obserr ! equivalence to handle character names equivalence(rstation_id,sidchr) !--- data statements data hdstr /'XOB YOB DHR'/ data oestr /'LOE'/ data qcstr /'LQM'/ data lunin / 13 / nreal=13 lob = obstype == 'goes_glm' ! . . . . ! Open, then read date from BUFR file open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) ! Initialization ntb = 0 nmsg = 0 disterrmax=-9999.0_r_kind allocate(cdata_all(nreal,maxobs),isort(maxobs)) isort = 0 cdata_all=zero nread=0 ntest=0 nvtest=0 nchanl=0 ilon=2 ilat=3 icntpnt=0 ! . . . . ! Big loop over glmbufr file : READING THE BUFR FILE loop_msg: do while (ireadmg(lunin,subset,idate)== 0) nmsg = nmsg+1 write(*,'(3a,i10)') 'subset=',subset,' cycle time =',idate loop_readsb: do while(ireadsb(lunin) == 0) ! use msg lookup table to decide which messages to skip ! use report id lookup table to only process matching reports ntb = ntb+1 ! Extract location and date information call ufbint(lunin,hdr,3,1,iret,hdstr) if(abs(hdr(2))>r90 .or. abs(hdr(1))>r360) cycle loop_readsb if(hdr(1) > r180)hdr(1)=hdr(1)-r360 if(hdr(1) < zero)hdr(1)=hdr(1)+r360 dlon_earth=hdr(1)*deg2rad dlat_earth=hdr(2)*deg2rad if (regional) then !-- WRF-ARW if (wrf_mass_regional) then call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) ! convert to rotated coordinates if(diagnostic_reg) then call txy2ll(dlon,dlat,rlon00,rlat00) ntest=ntest+1 cdist=sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* & (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00)) cdist=max(-one,min(cdist,one)) disterr=acos(cdist)*rad2deg disterrmax=max(disterrmax,disterr) end if if(outside) cycle loop_readsb ! check to see if outside regional domain endif ! wrf_mass_regional endif !if (regional) then ! Global if (.not. regional) then dlat = dlat_earth dlon = dlon_earth call grdcrd1(dlat,rlats,nlat,1) call grdcrd1(dlon,rlons,nlon,1) endif ! end global block if (offtime_data) then ! in time correction for observations to account for analysis ! time being different from obs file time. write(date,'( i10)') idate read (date,'(i4,3i2)') iy,im,idd,ihh idate5(1)=iy idate5(2)=im idate5(3)=idd idate5(4)=ihh idate5(5)=0 call w3fs21(idate5,minobs) ! obs ref time in seconds relative to historic date idate5(1)=iadate(1) idate5(2)=iadate(2) idate5(3)=iadate(3) idate5(4)=iadate(4) idate5(5)=0 call w3fs21(idate5,minan) ! analysis ref time in seconds relative to historic date ! Add obs reference time, then subtract analysis time to get obs time relative to analysis time_correction=float(minobs-minan)*r60inv else time_correction=zero end if timeobs=real(real(hdr(3),r_single),r_double) t4dv=timeobs + toff zeps=1.0e-8_r_kind if (t4dv -zeps) t4dv=zero if (t4dv>winlen.and.t4dvwinlen) cycle loop_readsb ! outside time window else if((real(abs(time)) > real(twindin)))cycle loop_readsb ! outside time window endif timex=time ! Note: An assessment of the GLM detection error is still a work in ! progress and at present, there are no clear metrics to assign a ! "quality mark" for an effective quality control procedure for ! lightning observations from this sensor. Nonetheless, the infrastructure ! for passing observations errors and quality control information to other ! routines within the GSI source code has been left in place. In addition, a ! temporarty undefined value (-.9999) has been assigned to the mnemonics ! corresponding to these variables in the lightning BUFR file. Developments ! for the calculation of observation errors and quality control (sanity checks) are ! expected in future upgrades to the "GOES/GLM variational lightning ! assimilation package." ! Extract observation error informatiom call ufbint(lunin,obserr,1,1,iret,oestr) loe=obserr(1,1) lmerr=loe ! Extract quality control information call ufbint(lunin,qcmark,1,1,iret,qcstr) lqm=qcmark(1,1) ! Data counter nread=nread+1 icntpnt=icntpnt+1 ndata=ndata+1 nodata=nodata+1 iout=ndata isort(icntpnt)=iout if (ndata > maxobs) then write(6,*)'READ_GOESGLM: ***WARNING*** ndata > maxobs for ',obstype ndata = maxobs end if ! Set usage variable usage = zero if (iuse_light(nlighttype) <= 0)usage=100._r_kind if (lob) then cdata_all(1,iout) =loe ! lightning observation error cdata_all(2,iout) =dlon ! grid relative longitude cdata_all(3,iout) =dlat ! grid relative latitude cdata_all(4,iout) =iout ! lightning obs cdata_all(5,iout) =rstation_id ! station id cdata_all(6,iout) =t4dv ! analysis time cdata_all(7,iout) =nlighttype ! type cdata_all(8,iout) =lmerr ! lightning max error cdata_all(9,iout) =lqm ! quality mark cdata_all(10,iout)=loe ! original lightning obs error loe cdata_all(11,iout)=usage ! usage parameter cdata_all(12,iout)=dlon_earth*rad2deg ! earth relative lon (degrees) cdata_all(13,iout)=dlat_earth*rad2deg ! earth relative lat (degrees) end if ! end loop on read line BUFR end do loop_readsb ! end of BUFR read loop enddo loop_msg !Uncomment if reading messeges in a loop ! . . . . ! Close unit to bufr file call closbf(lunin) ! 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,*) ' READ_GOESGLM: mix up in read_goesglm ,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 to the subroutine that transforms lightning strikes into lightning flash rate if(ndata /= 0) then !! count flash rate data and alocate temporary domain !! begin with the current number of strikes as the theoretical upper limit ndata_flash_h=ndata allocate(cdata_flash_h(nreal,ndata_flash_h)) call convert_to_flash_rate & (nreal,ndata,cdata_out,ndata_flash_h,cdata_flash_h,ndata_flash) deallocate(cdata_out) ndata=ndata_flash allocate(cdata_flash(nreal,ndata)) do i=1,ndata do k=1,nreal cdata_flash(k,i)=cdata_flash_h(k,i) end do end do deallocate(cdata_flash_h) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) cdata_flash !!!!!!!!! !!! Write lat, lon, time, lightning flash rate "superobs" into a file: !!! cdata_flash_h(2,iout),cdata_flash_h(3,iout),cdata_flash_h(4,iout), & !!! cdata_flash_h(6,iout)cdata_flash(4,:) deallocate(cdata_flash) else ! ndata=0 write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) cdata_out deallocate(cdata_out) end if !! if(ndata =/ 0) then 900 continue if(diagnostic_reg .and. ntest>0) write(6,*)'READ_GOESGLM: ',& 'ntest,disterrmax=',ntest,disterrmax if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_GOESGLM: ',& 'nvtest,vdisterrmax=',ntest,vdisterrmax if (ndata == 0) then write(6,*)'READ_GOESGLM: closbf(',lunin,') no data' endif close(lunin) ! End of routine return end subroutine read_goesglm !!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine convert_to_flash_rate & (nreal,ndata_strike,cdata_strike,ndata_flash_h,cdata_flash_h,ndata_flash) !$$$ documentation block ! . . . . ! subroutine: convert_to_flash_rate converts geo-located lightning strikes into ! lightning flash rate (#hits km-2 hr-1) ! prgmmr: zupanski ! org: CSU/CIRA, Data Assimilation group ! date: 2015-06-12 ! ! abstract: This subroutine does the following: !---- !- 1- counts the number of hits surrounding a GSI analysis grid point !- 2- calculates the flash rate averaged over time and area (hr*km**2) !- 3- finds the center of mass in terms of lon,lat and glon,glat !- 4- assigns a center of mass to be the flash rate observation point !---- ! ! program history log: ! 2015-07-01 apodaca - several updates in the calculation of lightning flashrate use kinds, only: r_single,r_kind,r_double,i_kind use constants, only: zero,one_tenth,one,deg2rad,& three,half,zero,& r10,r100,ten, r1000, rearth use gridmod, only: wrf_mass_regional,regional,nlon,nlat,& tll2xy,txy2ll use gridmod, only: lat2, lon2 use wrf_mass_guess_mod, only: ges_xlon, ges_xlat use gsi_4dvar, only: nhr_assimilation implicit none integer(i_kind),intent(in ) :: nreal,ndata_strike,ndata_flash_h integer(i_kind),intent(inout) :: ndata_flash real(r_kind),intent(inout),dimension(nreal,ndata_strike) :: cdata_strike real(r_kind),intent(inout),dimension(nreal,ndata_flash_h) :: cdata_flash_h real(r_kind),allocatable,dimension(:) :: gtim_central real(r_kind),allocatable,dimension(:) :: glon_central real(r_kind),allocatable,dimension(:) :: glat_central real(r_kind),allocatable,dimension(:) :: lon_central real(r_kind),allocatable,dimension(:) :: lat_central integer(i_kind),allocatable,dimension(:) :: ind_min integer(i_kind),allocatable,dimension(:) :: lcount real(r_kind) :: rearth2 real(r_kind) :: dtime,darea,cosine real(r_kind) :: darea_sum real(r_kind) :: delta_lon,delta_lat real(r_kind) :: lat_ref real(r_kind) :: xx,yy real(r_kind) :: dist2,dist_min integer(i_kind) :: ii0,jj0 integer(i_kind) :: ngridh integer(i_kind) :: index integer(i_kind) :: iobs,usage logical :: xflag,yflag real(r_kind) :: xbound,ybound integer(i_kind) :: nxdim,nydim integer(i_kind) :: icount ! Output files !---- !- unit grid lat-lon distances (Earth difference divided by grid difference) !- calculated from all averaged grid areas !!- note: this can be relaxed if the unit grid area (km) is known (i.e. darea) !! darea = (r*cos(lat)*dlon)*(dlat) rearth2=(rearth/r1000)**2 !! squared earth radius in km (need for darea calculation) if (ndata_strike>0) then darea_sum=0._r_kind do iobs=1,ndata_strike ii0=INT(cdata_strike(2,iobs)) jj0=INT(cdata_strike(3,iobs)) delta_lon=ges_xlon(jj0,ii0+1,1)-ges_xlon(jj0,ii0,1) delta_lat=ges_xlat(jj0+1,ii0,1)-ges_xlat(jj0,ii0,1) lat_ref =half*(cdata_strike(13,iobs)+cdata_strike(13,iobs-1)) cosine=cos(lat_ref*deg2rad) darea_sum=darea_sum+rearth2*cos(lat_ref*deg2rad)*(abs(delta_lon)*deg2rad)*& (abs(delta_lat)*deg2rad) end do !! do iobs=2,ndata_strike darea=darea_sum/float(ndata_strike) else !! ndata_strike=0 darea=zero end if !! if(ndata_strike>0) then dtime=float(nhr_assimilation) ! Regional if (regional) then !-- WRF-ARW if (wrf_mass_regional) then nxdim=lon2 nydim=lat2 endif ! wrf_mass_regional endif !if (regional) then ! Global if (.not. regional) then nxdim=nlon nydim=nlat endif ! end global block !!! Allocate new var for flash rate here (nxdim,nydim) ! asign zero to all points ! update the relevant points in the loop ngridh=nxdim*nydim allocate(gtim_central(1:ngridh)) allocate(glon_central(1:ngridh)) allocate(glat_central(1:ngridh)) allocate( lon_central(1:ngridh)) allocate( lat_central(1:ngridh)) allocate( lcount(1:ngridh)) lcount(:)=0 glon_central(:)=zero glat_central(:)=zero lon_central(:) =zero lat_central(:) =zero gtim_central(:)=zero do iobs=1,ndata_strike xx=cdata_strike(2,iobs) !! glon yy=cdata_strike(3,iobs) !! glat ii0=INT(cdata_strike(2,iobs)) jj0=INT(cdata_strike(3,iobs)) !! find lightning strikes near the (ii0,jj0) point xbound=float(ii0) ybound=float(jj0) xflag=(xx>xbound) .AND. (xxybound) .AND. (yy0) then glon_central(index)=glon_central(index)/float(lcount(index)) glat_central(index)=glat_central(index)/float(lcount(index)) lon_central(index)= lon_central(index)/float(lcount(index)) lat_central(index)= lat_central(index)/float(lcount(index)) endif !! if(lcount(index)>0) then enddo !! do index=1,ngridh !-- find the original index of the nearest strike (need for transfer of input obs) allocate(ind_min(1:ngridh)) ind_min(:)=-99 dist_min=1.e10_r_kind do iobs=1,ndata_strike xx=cdata_strike(2,iobs) !! glon yy=cdata_strike(3,iobs) !! glat ii0=INT(cdata_strike(2,iobs)) jj0=INT(cdata_strike(3,iobs)) index=(jj0-1)*nxdim+ii0 if (lcount(index)>0) then dist2=(xx-glon_central(index))**2+(yy-glat_central(index))**2 if (dist20) then enddo !! do iobs=1,ndata_strike !---- !---- Output !---- !-- count the non-zero flash rates and assign a temporary domain cdata_flash_h !-- Note: it is assumed that only non-zero flash rates are true observations icount=0 do index=1,ngridh if (lcount(index)>0) then icount=icount+1 cdata_flash_h( 1,icount)=cdata_strike( 1,ind_min(index)) cdata_flash_h( 2,icount)=glon_central(index) cdata_flash_h( 3,icount)=glat_central(index) if (darea>0._r_kind) then cdata_flash_h( 4,icount)=float(lcount(index))/(darea*dtime) else cdata_flash_h( 4,icount)=0. end if cdata_flash_h( 5,icount)=cdata_strike( 5,ind_min(index)) cdata_flash_h( 6,icount)=gtim_central(index) cdata_flash_h( 7,icount)=cdata_strike( 7,ind_min(index)) cdata_flash_h( 8,icount)=cdata_strike( 8,ind_min(index)) cdata_flash_h( 9,icount)=cdata_strike( 9,ind_min(index)) cdata_flash_h(10,icount)=cdata_strike(10,ind_min(index)) cdata_flash_h(11,icount)=usage cdata_flash_h(12,icount)=lon_central(index) cdata_flash_h(13,icount)=lat_central(index) endif !! if(lcount(index)>0) then enddo !! do index=1,ngridh ndata_flash=icount deallocate(ind_min) deallocate(lcount) deallocate(gtim_central) deallocate(glon_central) deallocate(glat_central) deallocate(lon_central) deallocate(lat_central) !----- ! End of routine return end subroutine convert_to_flash_rate !----- subroutine convert_time (date_old,date_new,nmax) !$$$ documentation block ! . . . . ! subroutine: convert_time ! prgmmr: k apodaca date: 2015-11-18 ! ! abstract: This subroutine performs a date conversion to deal with ! computer precission dependencies associated with a 10-digit ! float for the analysis date/time. !-- use kinds, only: r_kind,i_kind implicit none integer(i_kind), intent(in) :: nmax real(r_kind), intent(inout) :: date_old real(r_kind), intent(in) :: date_new integer(i_kind) :: i,sumidd integer(i_kind) :: idd,jdd,kdd real(r_kind), allocatable :: xdate(:) real(r_kind) :: dd,hh,ysumidd,xsumidd real(r_kind) :: xdd,xhh,ydate real(r_kind) :: xccyy allocate(xdate(1:nmax)) xdate(1:nmax-1) = date_old xdate(nmax) = date_new sumidd=0._r_kind do i=1,nmax xccyy = INT(1.0e-8_r_kind*xdate(i))*1.0e8_r_kind xdate(i) = INT(xdate(i))-xccyy jdd=INT(0.0001_r_kind*xdate(i)) idd=INT(xdate(i))-jdd*10000 ysumidd=float(idd) dd=float(INT(0.01_r_kind*ysumidd)) hh=ysumidd-dd*100._r_kind sumidd=sumidd+dd*24._r_kind+hh enddo !! do i=1,nmax xsumidd=float(sumidd)/nmax ysumidd=float(INT(xsumidd)) kdd=INT(xsumidd/24._r_kind) xdd=float(kdd) xhh=ysumidd-float(kdd)*24._r_kind ydate=float(jdd)*10000._r_kind+xdd*100._r_kind+xhh+xccyy date_old=ydate deallocate(xdate) end subroutine convert_time