subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,nobs) !$$$ subprogram documentation block ! . . . . ! subprogram: read_dbz read MRMS gridded QC'd radar reflectivity files in DART-like netcdf format ! ! abstract: Read and process MRMS gridded QC'd radar reflectivity (dBZ) ! observations in DART-like netcdf format. ! ! program history log: ! 2016-02-14 Y. Wang, Johnson, X. Wang - modify read_radar.f90 to read MRMS dbz in netcdf format ! in collaboration with Carley, POC: xuguang.wang@ou.edu ! 2019-02-27 D. Dowell : changed data_r_1d from real(r_kind) to real; added new array data_r_2d ! changed lon and lat to 2D arrays ! changed value for dbznoise ! ! program history log: ! ! input argument list: ! infile - file from which to read data ! lunout - unit to which to write data for further processing ! obstype - observation type to process ! ! output argument list: ! nread - number of radar reflectivity observations read ! ndata - number of radar reflectivity observations retained for further processing ! nodata - number of radar reflectivity observations retained for further processing ! sis - satellite/instrument/sensor indicator ! ! Variable Definitions: ! ! cdata_all - real - dim(maxdat,maxobs) - array holding all data for assimilation ! cstaid - char - radar station ide ! dlat - real - grid relative latitude of observation (grid units) ! dlon - real - grid relative longitude of observation (grid units) ! maxobs - int - max number of obs converted to no precip observations ! num_m2nopcp -int - number of missing obs ! num_missing - int - number of missing observations ! num_noise - int - number of rejected noise observations ! num_nopcp - int - number of noise obs converted to no precip observations ! numbadtime - int - number of elevations outside time window ! outside - logical - if observations are outside the domain -> true ! radartwindow - real - time window for radar observations (minutes) ! rmins_an - real - analysis time from reference date (minutes) ! rmins_ob - real - observation time from reference date (minutes) ! rstation_id - real - radar station id ! thisazimuthr - real - 90deg minues the actual azimuth and converted to radians ! thiserr - real - observation error ! thislat - real - latitude of observation, point ! thislon - real - longitude of observation, point ! thisrange - real - range of observation from radar ! thishgt - real - observation height, point ! this_stahgt - real - radar station height (meters about sea level) ! this_staid - char - radar station id ! thistiltr - real- radar tilt angle (radians) ! timeb - real - obs time (analyis relative minutes) ! dbzQC - real - reflectivity observation ! dbz_err - real - observation error of reflectivity ! height - real - height of observation ! lon - real - longitude of observation ! lat - real - latitude of observation ! utime - real - time for each observation point ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,r_double,i_kind,r_single use constants, only: zero,half,one,two,deg2rad,rad2deg, & one_tenth,r1000,r60,r60inv,r100,r400,grav_equator, & eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening use gridmod, only: tll2xy,nsig,nlat,nlon use obsmod, only: iadate,doradaroneob, & mintiltdbz,maxtiltdbz,minobrangedbz,maxobrangedbz,& static_gsi_nopcp_dbz,rmesh_dbz,zmesh_dbz use hybrid_ensemble_parameters,only : l_hyb_ens use obsmod,only: radar_no_thinning,missing_to_nopcp use convinfo, only: nconvtype,ctwind,icuse,ioctype use convthin, only: make3grids,map3grids,del3grids,use_all use jfunc, only: miter use mpimod, only: npe implicit none include 'netcdf.inc' ! Declare passed variables character(len=*),intent(in ) :: obstype,infile character(len=*),intent(in ) :: sis integer(i_kind) ,intent(in ) :: lunout integer(i_kind) ,intent(inout) :: nread,ndata,nodata real(r_kind),dimension(nlat,nlon,nsig),intent(in):: hgtl_full integer(i_kind),dimension(npe) ,intent(inout) :: nobs ! Declare local parameters real(r_kind),parameter:: r6 = 6.0_r_kind real(r_kind),parameter:: r360=360.0_r_kind integer(i_kind),parameter:: maxdat=18 character(len=4), parameter :: radid = 'XXXX' ! === Grid dbz data declaration real(r_single), allocatable, dimension(:) :: data_r_1d real(r_single), allocatable, dimension(:,:) :: data_r_2d real(r_kind), allocatable, dimension(:) :: height real(r_kind), allocatable, dimension(:,:) :: lon, lat real(r_single), allocatable, dimension(:,:,:) :: data_r_3d real(r_kind), allocatable, dimension(:,:,:) :: dbzQC integer(i_kind), parameter :: max_num_vars = 50, max_num_dims = 20 integer(i_kind) :: length, rcode, cdfid character( len = 20 ),dimension(max_num_vars) :: var_list integer(i_kind), dimension(max_num_vars) :: id_var, ndims integer(i_kind), dimension(max_num_dims) :: dimids, one_read integer(i_kind) :: natts, ivtype integer(i_kind) , dimension(max_num_vars, max_num_dims) :: dims logical :: if_input_exist integer(i_kind) :: ivar, var_num, sec70 !--Counters for diagnostics integer(i_kind) :: num_missing=0,num_nopcp=0, & !counts numbadtime=0, & num_m2nopcp=0, & num_noise=0,num_limmax=0 integer(i_kind) num_dbz2mindbz,imissing2nopcp integer(i_kind) :: ithin,zflag,nlevz,icntpnt,klon1,klat1,kk,klatp1,klonp1 real(r_kind) :: rmesh,xmesh,zmesh,dx,dy,dx1,dy1,w00,w01,w10,w11 real(r_kind), allocatable, dimension(:) :: zl_thin real(r_kind),dimension(nsig):: hges,zges real(r_kind) sin2,termg,termr,termrg,zobs,hgt integer(i_kind) ntmp,iout,iiout,ntdrvr_thin2 real(r_kind) crit1,timedif real(r_kind),parameter:: r16000 = 16000.0_r_kind logical :: luse integer(i_kind) maxout,maxdata integer(i_kind),allocatable,dimension(:):: isort !--General declarations integer(i_kind) :: ierror,i,j,k,nvol, & ikx,mins_an integer(i_kind) :: maxobs,nchanl,ilat,ilon,scount real(r_kind) :: thistiltr,thisrange,this_stahgt,thishgt real(r_kind) :: thisazimuthr,t4dv, & dlat,dlon,thiserr,thislon,thislat, & timeb real(r_kind) :: radartwindow real(r_kind) :: rmins_an real(r_kind),allocatable,dimension(:,:):: cdata_all real(r_double) rstation_id character(8) cstaid character(4) this_staid equivalence (this_staid,cstaid) equivalence (cstaid,rstation_id) logical :: outside real(r_kind) :: minobrange,maxobrange,mintilt,maxtilt real(r_kind) :: dbznoise=-10.0_r_kind ! dBZ obs must be >= dbznoise for assimilation logical :: l_limmax=.true. ! If true, observations > 60 dBZ are limited to be 60 dBZ. This is ! due to representativeness error associated with the model minobrange=minobrangedbz maxobrange=maxobrangedbz mintilt=mintiltdbz maxtilt=maxtiltdbz num_dbz2mindbz=0 write(6,*)'missing_to_nopcp is ',missing_to_nopcp write(6,*)'radar_no_thinning is ',radar_no_thinning !--------------------------------------------------------------------------------------! ! END OF ALL DECLARATIONS ! !--------------------------------------------------------------------------------------! write(6,*)'think in read_dbz static_gsi_nopcp_dbz is ', static_gsi_nopcp_dbz !-Check if reflectivity is in the convinfo file and extract necessary attributes ithin=1 !number of obs to keep per grid box if(radar_no_thinning) then ithin=-1 endif scount=0 ikx=0 do i=1,nconvtype if(trim(obstype) == trim(ioctype(i)) .and. abs(icuse(i))== 1) then ikx=i radartwindow=ctwind(ikx)*r60 !Time window units converted to minutes ! (default setting for dbz within convinfo is 0.05 hours) exit !Exit loop when finished with initial convinfo fields else if ( i==nconvtype ) then write(6,*) 'READ_dBZ: ERROR - OBSERVATION TYPE IS NOT PRESENT IN CONVINFO OR USE FLAG IS ZERO' write(6,*) 'READ_dBZ: ABORTTING read_dbz.f90 - NO REFLECTIVITY OBS READ!' return endif end do if (minobrange >= maxobrange) then write(6,*) 'MININMUM OB RANGE >= MAXIMUM OB RANGE FOR READING dBZ - PROGRAM STOPPING FROM READ_DBZ.F90' call stop2(400) end if !-next three values are dummy values for now nchanl=0 ilon=2 ilat=3 maxobs=50000000 !value taken from read_radar.f90 !--Allocate cdata_all array allocate(cdata_all(maxdat,maxobs),isort(maxobs)) rmesh=rmesh_dbz zmesh=zmesh_dbz maxout=0 maxdata=0 isort=0 ntdrvr_thin2=0 icntpnt=0 zflag=0 use_all=.true. if (ithin > 0) then write(6,*)'READ_RADAR_DBZ: ithin,rmesh :',ithin,rmesh use_all=.false. if(zflag == 0)then nlevz=nsig else nlevz=r16000/zmesh endif xmesh=rmesh call make3grids(xmesh,nlevz) ! call make3grids2(xmesh,nlevz) allocate(zl_thin(nlevz)) if (zflag == 1) then do k=1,nlevz zl_thin(k)=k*zmesh enddo endif write(6,*)'READ_RADAR_DBZ: xmesh, zflag, nlevz =', xmesh, zflag, nlevz endif !!end modified for thinning var_list(1:4) = (/ "reflectivity", "height ", "longitude ", "latitude "/) var_num = 4 print *, "read_Dbz.f90: open ",trim(infile) length = len_trim(infile) inquire(file=infile(1:length), exist=if_input_exist) fileopen: if (if_input_exist) then rcode = nf_open( infile(1:length), NF_NOWRITE, cdfid ) DO ivar = 1, var_num ! Check variable is in file, and get variable id: rcode = nf_inq_varid ( cdfid, trim(var_list(ivar)), id_var(ivar) ) if ( rcode /= 0 ) then write(6,FMT='(A,A)') & trim(var_list(ivar)), ' variable is not in input file' end if ! Get number of dimensions, and check of real type: dimids = 0 rcode = nf_inq_var( cdfid, id_var(ivar), trim(var_list(ivar)), ivtype, ndims(ivar), dimids, natts ) if ( ivtype /= 5 ) then write(6,FMT='(A,A)') trim(var_list(ivar)), ' variable is not real type' end if ! Get dimensions of field: dims(ivar,:) = 0 do i = 1, ndims(ivar) rcode = nf_inq_dimlen( cdfid, dimids(i), dims(ivar,i) ) end do END DO ! ivar allocate( dbzQC(dims(1,1),dims(1,2),dims(1,3)), height(dims(2,1)), & lon(dims(3,1),dims(3,2)), lat(dims(4,1),dims(4,2)) ) one_read = 1 do ivar = 1, var_num if( ivar == 1 )then allocate( data_r_3d(dims(ivar,1),dims(ivar,2),dims(ivar,3)) ) call ncvgt( cdfid, id_var(ivar), one_read, dims(ivar,:), data_r_3d, rcode ) dbzQC = data_r_3d else if( ivar == 2 )then allocate( data_r_1d(dims(ivar,1)) ) call ncvgt( cdfid, id_var(ivar), one_read, dims(ivar,:), data_r_1d, rcode ) height = data_r_1d deallocate( data_r_1d ) else allocate( data_r_2d(dims(ivar,1),dims(ivar,2)) ) call ncvgt( cdfid, id_var(ivar), one_read, dims(ivar,:), data_r_2d, rcode ) if( ivar == 3 )then lon = data_r_2d else if( ivar == 4 )then lat = data_r_2d end if deallocate( data_r_2d ) end if end do !-Obtain analysis time in minutes since reference date sec70 = 252460800 ! seconds since from 01/01/1970 call w3fs21(iadate,mins_an) !mins_an -integer number of mins snce 01/01/1978 rmins_an=mins_an !convert to real number ivar = 1 ILOOP : & do i = 1, dims(ivar,1) do j = 1, dims(ivar,2) do k = 1, dims(ivar,3) imissing2nopcp = 0 ! Missing data in the input file have the value -999.0 if( dbzQC(i,j,k) <= -900.0_r_kind ) then !--Extend no precip observations to missing data fields? ! May help suppress spurious convection if a problem. if (missing_to_nopcp .and. dbzQC(i,j,k) > -1000.0_r_kind ) then imissing2nopcp = 1 dbzQC(i,j,k) = static_gsi_nopcp_dbz num_m2nopcp = num_m2nopcp + 1 else num_missing = num_missing + 1 cycle end if end if if(miter /= 0 .and. (.not. l_hyb_ens) ) then ! For gsi 3DVar run if (l_limmax) then if ( dbzQC(i,j,k) > 60_r_kind ) then dbzQC(i,j,k) = 60_r_kind num_limmax = num_limmax + 1 end if end if end if if ( dbzQC(i,j,k) < static_gsi_nopcp_dbz ) then dbzQC(i,j,k) = static_gsi_nopcp_dbz num_dbz2mindbz = num_dbz2mindbz + 1 end if thishgt = height(k) ! unit : meter hgt = thishgt thislon = lon(i,j) thislat = lat(i,j) !-Check format of longitude and correct if necessary if(thislon>=r360) thislon=thislon-r360 if(thislon 0)then if(zflag == 0)then 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 hges(kk)=w00*hgtl_full(klat1 ,klon1 ,kk) + & w10*hgtl_full(klatp1,klon1 ,kk) + & w01*hgtl_full(klat1 ,klonp1,kk) + & w11*hgtl_full(klatp1,klonp1,kk) end do sin2 = sin(thislat)*sin(thislat) termg = grav_equator * & ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2)) termr = semi_major_axis /(one + flattening + grav_ratio - & two*flattening*sin2) termrg = (termg/grav)*termr do kk=1,nsig zges(kk) = (termr*hges(kk)) / (termrg-hges(kk)) zl_thin(kk)=zges(kk) end do endif zobs = hgt ntmp=ndata ! counting moved to map3gridS timedif=abs(t4dv) !don't know about this crit1 = timedif/r6+half call map3grids(1,zflag,zl_thin,nlevz,thislat,thislon,& zobs,crit1,ndata,iout,icntpnt,iiout,luse,.false.,.false.) maxout=max(maxout,iout) maxdata=max(maxdata,ndata) if (.not. luse) then ntdrvr_thin2=ntdrvr_thin2+1 cycle endif if(iiout > 0) isort(iiout)=0 if (ndata > ntmp) then nodata=nodata+1 endif isort(icntpnt)=iout else ndata =ndata+1 nodata=nodata+1 iout=ndata isort(icntpnt)=iout endif !!end modified for thinning thisazimuthr=0.0_r_kind this_staid=radid !Via equivalence in declaration, value is propagated ! to rstation_id used below. cdata_all(1,iout) = thiserr ! reflectivity obs error (dB) - inflated/adjusted cdata_all(2,iout) = dlon ! grid relative longitude cdata_all(3,iout) = dlat ! grid relative latitude cdata_all(4,iout) = thishgt ! obs absolute height (m) cdata_all(5,iout) = dbzQC(i,j,k) ! radar reflectivity factor cdata_all(6,iout) = thisazimuthr ! 90deg-azimuth angle (radians) cdata_all(7,iout) = timeb*r60inv ! obs time (analyis relative hour) cdata_all(8,iout) = ikx ! type cdata_all(9,iout) = thistiltr ! tilt angle (radians) cdata_all(10,iout)= this_stahgt ! station elevation (m) cdata_all(11,iout)= rstation_id ! station id cdata_all(12,iout)= icuse(ikx) ! usage parameter cdata_all(13,iout)= thislon*rad2deg ! earth relative longitude (degrees) cdata_all(14,iout)= thislat*rad2deg ! earth relative latitude (degrees) cdata_all(15,iout)= thisrange ! range from radar in m cdata_all(16,iout)= thiserr ! orginal error from convinfo file cdata_all(17,iout)= dbznoise ! noise threshold for reflectivity (dBZ) cdata_all(18,iout)= imissing2nopcp !=0, normal !=1, !values !converted !from !missing !values if(doradaroneob .and. (cdata_all(5,iout) > -99.0_r_kind) ) exit ILOOP end do ! k end do ! j end do ILOOP ! i if (.not. use_all) then deallocate(zl_thin) call del3grids endif !---all looping done now print diagnostic output write(6,*)'READ_dBZ: Reached eof on radar reflectivity file' write(6,*)'READ_dBZ: # volumes in input file =',nvol write(6,*)'READ_dBZ: # read in obs. number =',nread write(6,*)'READ_dBZ: # elevations outside time window =',numbadtime write(6,*)'READ_dBZ: # of noise obs to no precip obs =',num_nopcp write(6,*)'READ_dBZ: # of missing data to no precip obs =',num_m2nopcp write(6,*)'READ_dBZ: # of rejected noise obs =',num_noise write(6,*)'READ_dBZ: # of missing data =',num_missing write(6,*)'READ_dBZ: # changed to min dbz =',num_dbz2mindbz write(6,*)'READ_dBZ: # restricted to 60dBZ limit =',num_limmax !---Write observation to scratch file---! call count_obs(ndata,maxdat,ilat,ilon,cdata_all,nobs) write(lunout) obstype,sis,maxdat,nchanl,ilat,ilon write(lunout) ((cdata_all(k,i),k=1,maxdat),i=1,ndata) !---------------DEALLOCATE ARRAYS-------------! deallocate(cdata_all) else !fileopen write(6,*) 'READ_dBZ: ERROR OPENING RADAR REFLECTIVITY FILE: ',trim(infile),' IOSTAT ERROR: ',ierror, ' SKIPPING...' end if fileopen end subroutine read_dbz_nc