subroutine streamline_for_rjlist(cgrid) ! ! to compile: ! xlf90 streamline_for_rjlist.f -L/nwprod/lib -lw3_8 ! ! abstract: ! read in from original unformatted gsi diagnostic files ! and write out formatted streamlined files. Valid for ! gsi release 200609. implicit none integer(4),parameter::lun_t=10 integer(4),parameter::lun_q=11 integer(4),parameter::lun_ps=12 integer(4),parameter::lun_w=13 integer(4),parameter::lun_spd=14 character(60) cgrid character(8),allocatable,dimension(:):: cdiagbuf character(8),allocatable,dimension(:):: cprvstg character(8),allocatable,dimension(:):: csprvstg character(8) cstation character(3) otype integer(4) idate,nchar,nreal,i,ii,mype,lun real(4),allocatable,dimension(:,:)::rdiagbuf real(4) rlat,rlon,oberr,oberr2,ob,ob_model,ddiff, & dudiff,dvdiff,rmuse,xx,yy real(4) uob,uob_model,vob,vob_model,rfactor,qtflg,a0,b0 real(8) rlat8,rlon8,xx8,yy8 real(8) da8,alat18,elon18,elonv8,alatan8 integer(4) itype real (4) dtime logical fexist, r15887_diagfile_fmt namelist/diagfile_fmt/r15887_diagfile_fmt data r15887_diagfile_fmt/.true./ inquire(file='diagfile_fmt_input',exist=fexist) if (fexist) then open (55,file='diagfile_fmt_input',form='formatted') read(55,*) r15887_diagfile_fmt close(55) endif print*,'in streamline_for_rjlist: r15887_diagfile_fmt=',r15887_diagfile_fmt open(lun_t,file='stats_t_rj.out',form='formatted') !output file open(lun_q,file='stats_q_rj.out',form='formatted') !output file open(lun_ps,file='stats_ps_rj.out',form='formatted')!output file open(lun_w,file='stats_w_rj.out',form='formatted') !output file open(lun_spd,file='stats_spd_rj.out',form='formatted') !output file print*,'in streamline_for_rjlist: cgrid=',trim(cgrid) call proj_info(cgrid,da8,alat18,elon18,elonv8,alatan8) print*,'in streamline_for_rjlist: da8,alat18,elon18,& &elonv8,alatan8=',da8,alat18,elon18,elonv8,alatan8 open (7,file='diag_conv_for_rj.dat',form='unformatted') read(7,end=200) idate loop_read_obs: do read(7,end=200)otype,nchar,nreal,ii,mype ! print*,'in streamline_for_rjlist: idate=',idate ! print*,'in streamline_for_rjlist: otype,nchar,nreal,ii,mype=', & ! otype,nchar,nreal,ii,mype if (ii==0) then read(7,end=200) if (r15887_diagfile_fmt) read(7,end=200) cycle loop_read_obs endif allocate(cdiagbuf(ii)) allocate(cprvstg(ii)) allocate(csprvstg(ii)) allocate(rdiagbuf(nreal,ii)) if (r15887_diagfile_fmt) then read(7) cdiagbuf,rdiagbuf read(7) cprvstg,csprvstg else read(7) cdiagbuf,rdiagbuf,cprvstg,csprvstg endif if (otype(1:3)==' t') then lun=lun_t do i=1,ii cstation=cdiagbuf(i) itype=nint(rdiagbuf(1,i)) rlat=rdiagbuf(3,i) rlon=rdiagbuf(4,i) dtime=rdiagbuf(8,i) qtflg=rdiagbuf(10,i) if (rdiagbuf(11,i) .gt. 1.) then rmuse=rdiagbuf(11,i)*rdiagbuf(12,i) else rmuse=rdiagbuf(12,i) endif oberr=rdiagbuf(16,i) ob=rdiagbuf(17,i) ddiff=rdiagbuf(18,i) ob_model=rdiagbuf(17,i)-rdiagbuf(18,i) a0=0. b0=0. rlat8=rlat*1._8 rlon8=rlon*1._8 call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, & rlat8,rlon8,cgrid,xx8,yy8) xx=xx8 yy=yy8 oberr2=1.e10 if (oberr.gt.1.e-05) oberr2=1./oberr write(lun,'(a8)') cstation write(lun,*) itype,rlat,rlon,xx,yy,dtime,oberr2,ob,ob_model,a0,b0,rmuse enddo end if if (otype(1:3)==' q') then lun=lun_q do i=1,ii cstation=cdiagbuf(i) itype=nint(rdiagbuf(1,i)) rlat=rdiagbuf(3,i) rlon=rdiagbuf(4,i) dtime=rdiagbuf(8,i) qtflg=rdiagbuf(10,i) if (rdiagbuf(11,i) .gt. 1.) then rmuse=rdiagbuf(11,i)*rdiagbuf(12,i) else rmuse=rdiagbuf(12,i) endif oberr=rdiagbuf(16,i) ob=rdiagbuf(17,i) ddiff=rdiagbuf(18,i) ob_model=rdiagbuf(17,i)-rdiagbuf(18,i) a0=0. b0=0. rlat8=rlat*1._8 rlat8=rlat*1._8 rlon8=rlon*1._8 call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, & rlat8,rlon8,cgrid,xx8,yy8) xx=xx8 yy=yy8 oberr2=1.e10 if (oberr.gt.1.e-05) oberr2=1./oberr write(lun,'(a8)') cstation write(lun,*) itype,rlat,rlon,xx,yy,dtime,oberr2,ob,ob_model,a0,b0,rmuse enddo end if if (otype(2:3)=='ps') then lun=lun_ps do i=1,ii cstation=cdiagbuf(i) itype=nint(rdiagbuf(1,i)) rlat=rdiagbuf(3,i) rlon=rdiagbuf(4,i) dtime=rdiagbuf(8,i) qtflg=rdiagbuf(10,i) if (rdiagbuf(11,i) .gt. 1.) then rmuse=rdiagbuf(11,i)*rdiagbuf(12,i) else rmuse=rdiagbuf(12,i) endif oberr=rdiagbuf(16,i) ob=rdiagbuf(17,i) ddiff=rdiagbuf(18,i) ob_model=rdiagbuf(17,i)-rdiagbuf(18,i) a0=0. b0=0. rlat8=rlat*1._8 rlon8=rlon*1._8 call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, & rlat8,rlon8,cgrid,xx8,yy8) xx=xx8 yy=yy8 oberr2=1.e10 if (oberr.gt.1.e-05) oberr2=(1./oberr)*100. !convert to Pa ob=ob*100. !convert to Pa ob_model=ob_model*100.! convert to Pa write(lun,'(a8)') cstation write(lun,*) itype,rlat,rlon,xx,yy,dtime,oberr2,ob,ob_model,a0,b0,rmuse enddo end if if (otype(2:3)=='uv') then lun=lun_w do i=1,ii cstation=cdiagbuf(i) itype=nint(rdiagbuf(1,i)) rlat=rdiagbuf(3,i) rlon=rdiagbuf(4,i) dtime=rdiagbuf(8,i) qtflg=rdiagbuf(10,i) if (rdiagbuf(11,i) .gt. 1.) then rmuse=rdiagbuf(11,i)*rdiagbuf(12,i) else rmuse=rdiagbuf(12,i) endif oberr=rdiagbuf(16,i) uob=rdiagbuf(17,i) dudiff=rdiagbuf(18,i) uob_model=rdiagbuf(17,i)-rdiagbuf(18,i) vob=rdiagbuf(20,i) dvdiff=rdiagbuf(21,i) vob_model=rdiagbuf(20,i)-rdiagbuf(21,i) rfactor=rdiagbuf(23,i) rlat8=rlat*1._8 rlon8=rlon*1._8 call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, & rlat8,rlon8,cgrid,xx8,yy8) xx=xx8 yy=yy8 oberr2=1.e10 if (oberr.gt.1.e-05) oberr2=1./oberr write(lun,'(a8)') cstation write(lun,*) itype,rlat,rlon,xx,yy,dtime,oberr2,uob,uob_model,vob,vob_model,rmuse enddo end if if (otype(1:3)=='spd') then lun=lun_spd do i=1,ii cstation=cdiagbuf(i) itype=nint(rdiagbuf(1,i)) rlat=rdiagbuf(3,i) rlon=rdiagbuf(4,i) dtime=rdiagbuf(8,i) qtflg=rdiagbuf(10,i) if (rdiagbuf(11,i) .gt. 1.) then rmuse=rdiagbuf(11,i)*rdiagbuf(12,i) else rmuse=rdiagbuf(12,i) endif oberr=rdiagbuf(16,i) ob=rdiagbuf(17,i) ddiff=rdiagbuf(18,i) ob_model=rdiagbuf(17,i)-rdiagbuf(18,i) a0=0. b0=0. rfactor=rdiagbuf(20,i) rlat8=rlat*1._8 rlon8=rlon*1._8 call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, & rlat8,rlon8,cgrid,xx8,yy8) xx=xx8 yy=yy8 oberr2=1.e10 if (oberr.gt.1.e-05) oberr2=1./oberr write(lun,'(a8)') cstation write(lun,*) itype,rlat,rlon,xx,yy,dtime,oberr2,ob,ob_model,a0,b0,rmuse enddo end if deallocate(cdiagbuf,rdiagbuf) deallocate(cprvstg,csprvstg) enddo loop_read_obs 200 continue !123 format(a8,2x,i3,3(2x,f6.2),3(2x,f9.2)) !123 format(a8,2x,i3,2(2x,f6.2),4(2x,E11.4)) !123 format(a8,4x,i3,2x,2(2x,f6.2),3(2x,E11.4),2x,f5.1) !123 format(a8,4x,i3,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f5.1) 123 format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f5.1) close(lun_t) close(lun_q) close(lun_ps) close(lun_w) close(lun_spd) close(7) return end !*************************************************************** !note: !*************************************************************** ! if (otype(3:3)=='t') then !*************************************************************** ! rstation_id = data(id,i) ! cdiagbuf(ii) = station_id ! station id ! ! rdiagbuf(1,ii) = ictype(ikx) ! observation type ! rdiagbuf(2,ii) = icsubtype(ikx) ! observation subtype ! ! rdiagbuf(3,ii) = data(ilate,i) ! observation latitude (degrees) ! rdiagbuf(4,ii) = data(ilone,i) ! observation longitude (degrees) ! rdiagbuf(5,ii) = data(istnelv,i) ! station elevation (meters) ! rdiagbuf(6,ii) = prest ! observation pressure (hPa) ! rdiagbuf(7,ii) = data(iobshgt,i) ! observation height (meters) ! rdiagbuf(8,ii) = dtime ! obs time (hours relative to analysis time) ! rdiagbuf(9,ii) = data(iqc,i) ! input prepbufr qc or event mark ! rdiagbuf(10,ii) = data(iqt,i) ! setup qc or event mark (currently qtflg only) ! rdiagbuf(11,ii) = data(iuse,i) ! read_prepbufr data usage flag ! if(muse(i)) then ! rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) ! else ! rdiagbuf(12,ii) = -one ! endif ! ! err_input = data(ier2,i) ! err_adjst = data(ier,i) ! if (ratio_errors*error>tiny_r_kind) then ! err_final = one/(ratio_errors*error) ! else ! err_final = huge_single ! endif ! ! errinv_input = huge_single ! errinv_adjst = huge_single ! errinv_final = huge_single ! if (err_input>tiny_r_kind) errinv_input=one/err_input ! if (err_adjst>tiny_r_kind) errinv_adjst=one/err_adjst ! if (err_final>tiny_r_kind) errinv_final=one/err_final ! ! rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight ! rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (K**-1) ! rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error (K**-1) ! rdiagbuf(16,ii) = errinv_final ! final inverse observation error (K**-1) ! ! rdiagbuf(17,ii) = data(itob,i) ! temperature observation (K) ! rdiagbuf(18,ii) = ddiff ! obs-ges used in analysis (K) ! rdiagbuf(19,ii) = tob-tges ! obs-ges w/o bias correction (K) (future slot) !*************************************************************** !*************************************************************** ! if (otype(3:3)=='q') then !*************************************************************** ! rstation_id = data(id,i) ! cdiagbuf(ii) = station_id ! station id ! ! rdiagbuf(1,ii) = ictype(ikx) ! observation type ! rdiagbuf(2,ii) = icsubtype(ikx) ! observation subtype ! ! rdiagbuf(3,ii) = data(ilate,i) ! observation latitude (degrees) ! rdiagbuf(4,ii) = data(ilone,i) ! observation longitude (degrees) ! rdiagbuf(5,ii) = data(istnelv,i) ! station elevation (meters) ! rdiagbuf(6,ii) = presq ! observation pressure (hPa) ! rdiagbuf(7,ii) = data(iobshgt,i) ! observation height (meters) ! rdiagbuf(8,ii) = dtime ! obs time (hours relative to analysis time) ! ! rdiagbuf(9,ii) = data(iqc,i) ! input prepbufr qc or event mark ! rdiagbuf(10,ii) = rmiss_single ! setup qc or event mark ! rdiagbuf(11,ii) = data(iuse,i) ! read_prepbufr data usage flag ! if(muse(i)) then ! rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) ! else ! rdiagbuf(12,ii) = -one ! endif ! ! err_input = data(ier2,i)*qsges ! convert rh to q ! err_adjst = data(ier,i)*qsges ! convert rh to q ! if (ratio_errors*error>tiny_r_kind) then ! err_final = one/(ratio_errors*error) ! else ! err_final = huge_single ! endif ! ! errinv_input = huge_single ! errinv_adjst = huge_single ! errinv_final = huge_single ! if (err_input>tiny_r_kind) errinv_input = one/err_input ! if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst ! if (err_final>tiny_r_kind) errinv_final = one/err_final ! ! rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight ! rdiagbuf(14,ii) = errinv_input ! prepbufr inverse observation error ! rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error ! rdiagbuf(16,ii) = errinv_final ! final inverse observation error ! ! rdiagbuf(17,ii) = data(iqob,i) ! observation ! rdiagbuf(18,ii) = ddiff ! obs-ges used in analysis ! rdiagbuf(19,ii) = qob-qges ! obs-ges w/o bias correction (future slot) ! ! rdiagbuf(20,ii) = qsges ! guess saturation specific humidity ! !*************************************************************** !*************************************************************** ! if (otype(2:3)=='ps') then !*************************************************************** ! rstation_id = data(id,i) ! cdiagbuf(ii) = station_id ! station id ! ! rdiagbuf(1,ii) = ictype(ikx) ! observation type ! rdiagbuf(2,ii) = icsubtype(ikx) ! observation subtype ! rdiagbuf(3,ii) = data(ilate,i) ! observation latitude (degrees) ! rdiagbuf(4,ii) = data(ilone,i) ! observation longitude (degrees) ! rdiagbuf(5,ii) = data(istnelv,i) ! station elevation (meters) ! rdiagbuf(6,ii) = data(ipres,i)*r10 ! observation pressure (hPa) ! rdiagbuf(7,ii) = dhgt ! observation height (meters) ! rdiagbuf(8,ii) = dtime ! obs time (hours relative to analysis time) ! ! rdiagbuf(9,ii) = data(iqc,i) ! input prepbufr qc or event mark ! rdiagbuf(10,ii) = rmiss_single ! setup qc or event mark ! rdiagbuf(11,ii) = data(iuse,i) ! read_prepbufr data usage flag ! if(muse(i)) then ! rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) ! else ! rdiagbuf(12,ii) = -one ! endif ! ! pob = pob*r10 ! pges = pges*r10 ! pgesorig = pgesorig*r10 ! ! err_input = data(ier2,i)*r10 ! r10 converts cb to mb ! err_adjst = data(ier,i)*r10 ! if (ratio_errors*error/r10>tiny_r_kind) then ! err_final = r10/(ratio_errors*error) ! else ! err_final = huge_single ! endif ! ! errinv_input = huge_single ! errinv_adjst = huge_single ! errinv_final = huge_single ! if (err_input>tiny_r_kind) errinv_input = one/err_input ! if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst ! if (err_final>tiny_r_kind) errinv_final = one/err_final ! ! rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight ! rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (hPa**-1) ! rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error (hPa**-1) ! rdiagbuf(16,ii) = errinv_final ! final inverse observation error (hPa**-1) ! ! rdiagbuf(17,ii) = pob ! surface pressure observation (hPa) ! rdiagbuf(18,ii) = pob-pges ! obs-ges used in analysis (coverted to hPa) ! rdiagbuf(19,ii) = pob-pgesorig ! obs-ges w/o adjustment to guess surface pressure (hPa) !*************************************************************** !*************************************************************** ! if (otype(2:3)=='uv') then !*************************************************************** ! rstation_id = data(id,i) ! cdiagbuf(ii) = station_id ! station id ! ! rdiagbuf(1,ii) = ictype(ikx) ! observation type ! rdiagbuf(2,ii) = icsubtype(ikx) ! observation subtype ! ! rdiagbuf(3,ii) = data(ilate,i) ! observation latitude (degrees) ! rdiagbuf(4,ii) = data(ilone,i) ! observation longitude (degrees) ! rdiagbuf(5,ii) = data(ielev,i) ! station elevation (meters) ! rdiagbuf(6,ii) = presw ! observation pressure (hPa) ! rdiagbuf(7,ii) = data(ihgt,i) ! observation height (meters) ! rdiagbuf(8,ii) = dtime ! obs time (hours relative to analysis time) ! ! rdiagbuf(9,ii) = data(iqc,i) ! input prepbufr qc or event mark ! rdiagbuf(10,ii) = rmiss_single ! setup qc or event mark ! rdiagbuf(11,ii) = data(iuse,i) ! read_prepbufr data usage flag ! if(muse(i)) then ! rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) ! else ! rdiagbuf(12,ii) = -one ! endif ! ! err_input = data(ier2,i) ! err_adjst = data(ier,i) ! if (ratio_errors*error>tiny_r_kind) then ! err_final = one/(ratio_errors*error) ! else ! err_final = huge_single ! endif ! ! errinv_input = huge_single ! errinv_adjst = huge_single ! errinv_final = huge_single ! if (err_input>tiny_r_kind) errinv_input = one/err_input ! if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst ! if (err_final>tiny_r_kind) errinv_final = one/err_final ! ! rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight ! rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (m/s)**-1 ! rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error (m/s)**-1 ! rdiagbuf(16,ii) = errinv_final ! final inverse observation error (m/s)**-1 ! ! rdiagbuf(17,ii) = data(iuob,i) ! u wind component observation (m/s) ! rdiagbuf(18,ii) = dudiff ! u obs-ges used in analysis (m/s) ! rdiagbuf(19,ii) = uob-ugesin ! u obs-ges w/o bias correction (m/s) (future slot) ! ! rdiagbuf(20,ii) = data(ivob,i) ! v wind component observation (m/s) ! rdiagbuf(21,ii) = dvdiff ! v obs-ges used in analysis (m/s) ! rdiagbuf(22,ii) = vob-vgesin ! v obs-ges w/o bias correction (m/s) (future slot) ! ! rdiagbuf(23,ii) = factw ! 10m wind reduction factor !*************************************************************** !*************************************************************** ! if (otype(1:3)=='spd') then !*************************************************************** ! rstation_id = data(id,i) ! cdiagbuf(ii) = station_id ! station id ! ! rdiagbuf(1,ii) = ictype(ikx) ! observation type ! rdiagbuf(2,ii) = icsubtype(ikx) ! observation subtype ! ! rdiagbuf(3,ii) = data(ilate,i) ! observation latitude (degrees) ! rdiagbuf(4,ii) = data(ilone,i) ! observation longitude (degrees) ! rdiagbuf(5,ii) = data(istnelv,i) ! station elevation (meters) ! rdiagbuf(6,ii) = presw ! observation pressure (hPa) ! rdiagbuf(7,ii) = data(ihgt,i) ! observation height (meters) ! rdiagbuf(8,ii) = dtime ! obs time (hours relative to analysis time) ! ! rdiagbuf(9,ii) = data(iqc,i) ! input prepbufr qc or event mark ! rdiagbuf(10,ii) = rmiss_single ! setup qc or event mark ! rdiagbuf(11,ii) = data(iuse,i) ! read_prepbufr data usage flag ! if(muse(i)) then ! rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) ! else ! rdiagbuf(12,ii) = -one ! endif ! ! spdob0 = sqrt(data(iuob,i)*data(iuob,i)+data(ivob,i)*data(ivob,i)) ! err_input = data(ier2,i) ! err_adjst = data(ier,i) ! if (ratio_errors*error>tiny_r_kind) then ! err_final = one/(ratio_errors*error) ! else ! err_final = huge_single ! endif ! ! errinv_input = huge_single ! errinv_adjst = huge_single ! errinv_final = huge_single ! if (err_input>tiny_r_kind) errinv_input = one/err_input ! if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst ! if (err_final>tiny_r_kind) errinv_final = one/err_final ! ! rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight ! rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (m/s)**-1 ! rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error (m/s)**-1 ! rdiagbuf(16,ii) = errinv_final ! final inverse observation error (m/s)**-1 ! ! rdiagbuf(17,ii) = spdob ! wind speed observation (m/s) ! rdiagbuf(18,ii) = ddiff ! obs-ges used in analysis (m/s) ! rdiagbuf(19,ii) = spdob0-spdges ! obs-ges w/o bias correction (m/s) (future slot) ! ! rdiagbuf(20,ii) = factw ! 10m wind reduction factor !***************************************************************