diff --git a/cmake/Modules/FindBUFR.cmake b/cmake/Modules/FindBUFR.cmake index 58527743bc..a2d70b92bd 100644 --- a/cmake/Modules/FindBUFR.cmake +++ b/cmake/Modules/FindBUFR.cmake @@ -11,8 +11,8 @@ endif() set( NO_DEFAULT_PATH ) if(NOT BUILD_BUFR ) - if(DEFINED ENV{BUFR_LIBd} ) - set(BUFR_LIBRARY $ENV{BUFR_LIBd} ) + if(DEFINED ENV{BUFR_LIBd_DA} ) + set(BUFR_LIBRARY $ENV{BUFR_LIBd_DA} ) message("BUFR library ${BUFR_LIBRARY} set via Environment variable") else() find_library( BUFR_LIBRARY diff --git a/cmake/Modules/FindCORELIBS.cmake b/cmake/Modules/FindCORELIBS.cmake index 711864ebce..9aab567113 100644 --- a/cmake/Modules/FindCORELIBS.cmake +++ b/cmake/Modules/FindCORELIBS.cmake @@ -102,8 +102,8 @@ else() set( w3nco "w3nco${libsuffix}") endif() if(NOT BUILD_BUFR ) - if(DEFINED ENV{BUFR_LIBd} ) - set(BUFR_LIBRARY $ENV{BUFR_LIBd} ) + if(DEFINED ENV{BUFR_LIBd_DA} ) + set(BUFR_LIBRARY $ENV{BUFR_LIBd_DA} ) else() find_library( BUFR_LIBRARY NAMES libbufr.a libbufr_d_64.a libbufr_i4r8.a libbufr_v${BUFR_VER}_d_64.a diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90 index 4b98d904b6..a7f1ab0e4d 100644 --- a/src/gsi/cplr_gfs_ensmod.f90 +++ b/src/gsi/cplr_gfs_ensmod.f90 @@ -822,6 +822,10 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi if (k3u==0.or.k3v==0.or.k3t==0.or.k3q==0.or.k3oz==0) & write(6,'(" WARNING, problem with one of k3-")') + do k=1,nsig + call fillpoles_sv_(temp3(:,:,k,k3u),temp3(:,:,k,k3v),nlon,nlat,clons,slons) + end do + ! convert T to Tv: postpone this calculation ! temp3(:,:,:,k3t)=temp3(:,:,:,k3t)*(one+fv*temp3(:,:,:,k3q)) diff --git a/src/gsi/genqsat.f90 b/src/gsi/genqsat.f90 index 37927f1847..2719ed28f8 100644 --- a/src/gsi/genqsat.f90 +++ b/src/gsi/genqsat.f90 @@ -53,7 +53,7 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative) !$$$ use kinds, only: r_kind,i_kind use constants, only: xai,tmix,xb,omeps,eps,xbi,one,zero,& - xa,psat,ttp,half,one_tenth + xa,psat,ttp,half,one_tenth,qmin use derivsmod, only: qgues,dqdt,dqdrh,dqdp use jfunc, only: pseudo_q2 use gridmod, only: wrf_nmm_regional,wrf_mass_regional,nems_nmmb_regional,aeta2_ll,regional,cmaq_regional @@ -161,6 +161,7 @@ subroutine genqsat(qsat,tsen,prsl,lat2,lon2,nsig,ice,iderivative) end if es2=min(es,esmax) qsat(i,j,k) = eps * es2 / (pw - omeps * es2) + qsat(i,j,k) = max(qmin,qsat(i,j,k)) if(iderivative > 0)then if(es <= esmax .and. iderivative == 2)then diff --git a/src/gsi/getcount_bufr.f90 b/src/gsi/getcount_bufr.f90 index 1af117dfc2..40b68cf85c 100644 --- a/src/gsi/getcount_bufr.f90 +++ b/src/gsi/getcount_bufr.f90 @@ -33,7 +33,7 @@ subroutine getcount_bufr(inpfile,nmsg,nsub) lunit=get_lun() nsub=0;nmsg=0 - call closbf(lunit) +! call closbf(lunit) open(lunit,file=trim(inpfile),form='unformatted') call openbf(lunit,'IN',lunit) do while(ireadmg(lunit,subset,idate) >=0) diff --git a/src/gsi/glbsoi.f90 b/src/gsi/glbsoi.f90 index 222b49e832..d9dbee7b3b 100644 --- a/src/gsi/glbsoi.f90 +++ b/src/gsi/glbsoi.f90 @@ -169,7 +169,7 @@ subroutine glbsoi ! Declare local variables logical laltmin - integer(i_kind) jiterlast + integer(i_kind) jiterlast,lunix,lunit real(r_kind) :: zgg,zxy character(len=12) :: clfile logical print_verbose @@ -210,8 +210,17 @@ subroutine glbsoi ! Set cost function call create_jfunc + lunit=0 + lunix=0 + call isetprm('MXMSGL',400000) + call isetprm('MAXSS',250000) +! Initialize bufr read on all processors (so that exitbufr works, lunit and +! lunix are dummys) + call openbf(lunit,'FIRST',lunix) ! Read observations and scatter call observer_set +! release bufr memory + call exitbufr ! cloud analysis if(i_gsdcldanal_type==6 .or. i_gsdcldanal_type==3) then diff --git a/src/gsi/hdraobmod.f90 b/src/gsi/hdraobmod.f90 new file mode 100644 index 0000000000..d7dc6a810e --- /dev/null +++ b/src/gsi/hdraobmod.f90 @@ -0,0 +1,1315 @@ +module hdraobmod +!$$$ module documentation block +! . . . . +! module: obsmod +! prgmmr: derber org: np23 date: 2003-09-25 +! +! abstract: This module contains variables and arrays pertinent for +! observational data. +! +! program history log: +! 2021-03-18 derber +! +! def nhdps - number of high resolution surface pressure stations +! def nhduv - number of high resolution uv stations +! def nhdt - number of high resolution temperature stations +! def nhdq - number of high resolution moisture stations +! def hdpslist - list of high resolution surface pressure stations +! def hduvlist - list of high resolution uv stations +! def hdtlist - list of high resolution temperature stations +! def hdqlist - list of high resolution moisture stations +! +! attributes: +! langauge: f90 +! machine: +! +!$$$ end documentation block + + use kinds, only: i_kind + implicit none + +! set default as private + private +! set subroutines and functions to public + public :: read_hdraob + public :: nhdt,nhdq,nhduv,nhdps + public :: nodet,nodeq,nodeps,nodeuv + public :: hdtlist,hdqlist,hduvlist,hdpslist + + integer(i_kind) :: nhdt,nhdq,nhduv,nhdps + integer(i_kind) :: nodet,nodeq,nodeps,nodeuv + integer(i_kind),allocatable,dimension(:) :: hdpslist,hduvlist,hdtlist,hdqlist + +contains + subroutine read_hdraob(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& + prsl_full,hgtl_full,nobs,nrec_start) +!$$$ subprogram documentation block +! . . . . +! subprogram: read_hdraob read obs from hdraob file +! prgmmr: derber org: np22 date: 2020-02-12 +! +! abstract: This routine reads high resolution raob data found in the hdraob +! file. Specific observation types read by this routine +! include surface pressure, temperature, winds (components +! and speeds), and moisture. +! +! When running the gsi in regional mode, the code only +! retains those observations that fall within the regional +! domain +! +! program history log: +! 2020-02-24 derber + +! 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 +! prsl_full- 3d pressure on full domain grid +! hgtl_full- 3d height on full domain grid +! nrec_start - number of subsets without useful information +! +! 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 +! nobs - array of observations on each subdomain for each processor +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_single,r_kind,r_double,i_kind + use constants, only: zero,one_tenth,one,deg2rad,half,& + rad2deg,tiny_r_kind,huge_r_kind,& + r60inv,r10,r100,r2000,eps,omeps + use constants,only: grav + use gridmod, only: diagnostic_reg,regional,nlon,nlat,nsig,& + tll2xy,rotate_wind_ll2xy,rotate_wind_xy2ll,& + rlats,rlons,twodvar_regional,fv3_regional + use convinfo, only: nconvtype,ctwind, & + ncmiter,ncgroup,ncnumgrp,icuse,ictype,icsubtype,ioctype + + use obsmod, only: iadate,oberrflg,perturb_obs,perturb_fact,ran01dom + use obsmod, only: blacklst,offtime_data + use converr,only: etabl + use converr_ps,only: etabl_ps,isuble_ps,maxsub_ps + use converr_q,only: etabl_q,isuble_q,maxsub_q + use converr_t,only: etabl_t,isuble_t,maxsub_t + use converr_uv,only: etabl_uv,isuble_uv,maxsub_uv + use convb_ps,only: btabl_ps + use convb_q,only: btabl_q + use convb_t,only: btabl_t + use convb_uv,only: btabl_uv + use gsi_4dvar, only: time_4dvar,winlen + use qcmod, only: errormod,errormod_hdraob,njqc + use mpimod, only: mype +! use blacklist, only : blacklist_read,blacklist_destroy + use ndfdgrids,only: adjust_error + use deter_sfc_mod, only: deter_sfc2 + use mpimod, only: npe + use gsi_io, only: verbose + + implicit none + +! Declare passed variables + character(len=*) ,intent(in ) :: infile,obstype + character(len=20) ,intent(in ) :: sis + integer(i_kind) ,intent(in ) :: lunout,nrec_start + integer(i_kind) ,intent(inout) :: nread,ndata,nodata + integer(i_kind),dimension(npe) ,intent(inout) :: nobs + real(r_kind) ,intent(in ) :: twindin + real(r_kind),dimension(nlat,nlon,nsig),intent(in ) :: prsl_full,hgtl_full + +! Declare local parameters + real(r_kind),parameter:: r1_2 = 1.2_r_kind + real(r_kind),parameter:: r50 = 50.0_r_kind + real(r_kind),parameter:: r90 = 90.0_r_kind + real(r_kind),parameter:: r360 = 360.0_r_kind + +! integer(i_kind),parameter:: mxtb=5000000 +! integer(i_kind),parameter:: nmsgmax=100000 ! max message count + integer(i_kind), parameter:: maxlevs=20000 + +! Declare local variables + logical tob,qob,uvob,psob + logical outside + logical,allocatable,dimension(:,:):: lmsg ! set true when convinfo entry id found in a message + + character(80) hdstr,hdstr2,hdstr3,levstr,hdtypestr + character(80) obstr + character(10) date + character(8) subset + character(8) prvstr,sprvstr + character(8) c_station_id,id_station + character(8) stnid + character(8) stntype + + integer(i_kind) ireadmg,ireadsb + integer(i_kind) lunin,i,maxobs,j,idomsfc,nmsgmax,mxtb + integer(i_kind) kk,klon1,klat1,klonp1,klatp1 + integer(i_kind) nc,nx,ntread,itx,ii,ncsave + integer(i_kind) ihh,idd,idate,iret,im,iy,k,levs + integer(i_kind) kx,nreal,nchanl,ilat,ilon + integer(i_kind) lim_qm + integer(i_kind) iout + integer(i_kind) irec + integer(i_kind) ntest,nvtest + integer(i_kind) kl,k1,k2,k1_ps,k1_q,k1_t,k1_uv,k2_q,k2_t,k2_uv,k2_ps + integer(i_kind) itypex,id + integer(i_kind) minobs,minan + integer(i_kind) ntb,ntmatch,ncx,cat + integer(i_kind) nmsg ! message index + integer(i_kind) jj,start,next,ncount_ps,ncount_q,ncount_uv,ncount_t + integer(i_kind),dimension(5):: idate5 + integer(i_kind),dimension(nconvtype)::ntxall + integer(i_kind),dimension(nconvtype+1)::ntx + integer(i_kind),allocatable,dimension(:,:):: tab + integer(i_kind) igroup,istation + integer(i_kind) ierr_ps,ierr_q,ierr_t,ierr_uv ! the position of error table collum + integer(i_kind),dimension(maxlevs):: pqm,qqm,tqm,wqm + integer(i_kind) nstations + integer(i_kind),allocatable,dimension(:):: list_stations + real(r_kind) time,timeobs,toff,t4dv,zeps + real(r_kind) qtflg,ediff,usage,ediff_ps,ediff_q,ediff_t,ediff_uv + real(r_kind) u0,v0,uob,vob,dx,dy,dx1,dy1,w00,w10,w01,w11 + real(r_kind) qoe,qobcon,dlnpob,ppb,poe,obheight + real(r_kind) toe,woe,errout,oelev,dlat,dlon,dlat_earth,dlon_earth + real(r_kind) dlat_earth_deg,dlon_earth_deg + real(r_kind) elev,stnelev + real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00 + real(r_kind) vdisterrmax,u00,v00 + real(r_kind) rminobs,rminan,es,dummy + real(r_kind) del,terrmin,werrmin,perrmin,qerrmin,del_ps,del_q,del_t,del_uv + real(r_kind) tsavg,ff10,sfcr,zz,deltat + real(r_kind) time_correction,fact + real(r_kind),dimension(nsig):: presl,hgtl,lnpresl + real(r_kind),dimension(nsig-1):: dpres + real(r_kind),dimension(maxlevs)::plevs + real(r_kind),allocatable,dimension(:,:):: cdata_all,cdata_out + real(r_kind) :: zob,tref + + real(r_double) rstation_id,r_station + real(r_double),dimension(2):: hdr + real(r_double),dimension(9):: hdr2 + real(r_double),dimension(1):: hdr3 + real(r_double),dimension(1):: hdrtype + real(r_double),dimension(2,maxlevs):: levdat + real(r_double),dimension(8,maxlevs):: var_jb,obserr + real(r_double),dimension(8,maxlevs):: obsdat + + logical newstation,toocold + + +! equivalence to handle character names + equivalence(rstation_id,c_station_id) + equivalence(r_station,id_station) + +! data statements + data hdstr /'WMOB WMOS'/ + data hdstr2 /'YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH' / + data hdtypestr /'BUHD' / + data hdstr3 /'HEIT'/ + data obstr /'LTDS LATDH LONDH GP10 WSPD WDIR TMDB TMDP'/ + data levstr /'PRLC GP07'/ + data lunin / 13 / + !* for match loction station and time +! character(7*2000) cstn_idtime,cstn_idtime2 +! character(7) stn_idtime(2000),stn_idtime2(2000) +! equivalence (stn_idtime(1),cstn_idtime) +! equivalence (stn_idtime2(1),cstn_idtime2) +! integer :: ii1,atmp,btmp,mytimeyy,ibyte +! character(4) stid +! real(8) :: rval +! character(len=8) :: cval +! equivalence (rval,cval) +! character(7) flnm + + logical print_verbose,descend + + print_verbose=.false. + if(verbose) print_verbose=.true. + +! Initialize variables + + tob = obstype == 't' + uvob = obstype == 'uv' + qob = obstype == 'q' + psob = obstype == 'ps' + + nstations=0 + nreal=0 + if(tob)then + nreal=25 + else if(uvob) then + nreal=26 + else if(psob) then + nreal=20 + else if(qob) then + nreal=26 + else + write(6,*) ' illegal obs type in READ_HDRAOB ',obstype + call stop2(94) + end if + + if(perturb_obs .and. (tob .or. psob .or. qob))nreal=nreal+1 + if(perturb_obs .and. uvob )nreal=nreal+2 + + +! if (blacklst) call blacklist_read(obstype) + + lim_qm=3 + terrmin=half + werrmin=one + perrmin=0.3_r_kind + qerrmin=0.05_r_kind +!------------------------------------------------------------------------ + ntread=1 + ntmatch=0 + ntx(ntread)=0 + ntxall=0 + var_jb=zero + + + do nc=1,nconvtype + if(trim(ioctype(nc)) == trim(obstype))then + ntmatch=ntmatch+1 + ntxall(ntmatch)=nc + end if + end do + + irec = 0 +!! get message and subset counts + + call getcount_bufr(infile,nmsgmax,mxtb) + allocate(lmsg(nmsgmax,ntread),tab(mxtb,3)) + + lmsg = .false. + maxobs=0 + tab=0 + nmsg=0 + ntb = 0 + ncount_ps=0;ncount_q=0;ncount_t=0;ncount_uv=0 + +! Open, then read date from bufr data +! call closbf(lunin) + open(lunin,file=trim(infile),form='unformatted') + call openbf(lunin,'IN',lunin) + call datelen(10) + +! This is the message creating error message + msg_report: do while (ireadmg(lunin,subset,idate) == 0) + irec = irec + 1 + if(irec < nrec_start) cycle msg_report + +! Time offset + if(nmsg == 0) call time_4dvar(idate,toff) + nmsg=nmsg+1 + if (nmsg>nmsgmax) then + write(6,*)'READ_HDRAOB: messages exceed maximum ',nmsgmax + call stop2(50) + endif + loop_report: do while (ireadsb(lunin) == 0) + ntb = ntb+1 + if (ntb>mxtb) then + write(6,*)'READ_HDRAOB: reports exceed maximum ',mxtb + call stop2(50) + endif + + +! call ufbint(lunin,hdr,2,1,iret,hdstr) +! igroup=hdr(1) +! istation=hdr(2) +! write(6,*) igroup,istation +! call ufbint(lunin,obsdat,13,maxlevs,levs,obstr) +! Check for blacklisting of station ID +! if (blacklst .and. ibcnt > 0) then +! stnid = transfer(hdr(4),stnid) +! do i = 1,ibcnt +! if( kx == blkkx(i) .and. stnid == blkstns(i) ) then +! write(6,*)'READ_HDRAOB: blacklist station ',stnid, & +! 'for obstype ',trim(obstype),' and kx=',kx +! cycle loop_report +! endif +! enddo +! endif + + call ufbint(lunin,levdat,2,maxlevs,levs,levstr) + if(levs > maxlevs)write(6,*) ' not enough levels increase maxlevs ',levs,maxlevs + call ufbint(lunin,obsdat,8,maxlevs,levs,obstr) + deltat=0. + descend=.false. + if(obsdat(1,1) < 1.e6 .and. obsdat(1,levs) < 1.e6)then + deltat=obsdat(1,levs)-obsdat(1,1) + if(deltat < -1.) descend=.true. + end if +! Extract type information + if(psob .or. tob .or. qob)then + if(.not. descend) then + kx=119 + else + kx=118 + if(psob) cycle loop_report + end if + else if(uvob)then + if(.not. descend) then + if(levdat(1,2) > 1.e8)then + kx=218 + else + kx=219 + end if + else + if(levdat(1,2) > 1.e8)then + cycle loop_report ! should not have descent obs with no pressure + else + kx=217 + end if + end if + end if +! Match ob to proper convinfo type + ncsave=0 + matchloop:do ncx=1,ntmatch + nc=ntxall(ncx) + if (kx == ictype(nc))then + + ncsave=nc + exit matchloop + end if + + end do matchloop + +! write(6,*) ' levs =',levs + maxobs=maxobs+max(1,levs) + nx=ictype(ncsave) + tab(ntb,1)=ncsave + tab(ntb,2)=nx + tab(ntb,3)=levs + + end do loop_report + enddo msg_report + if (nmsg==0) then + call closbf(lunin) + close(lunin) + if(print_verbose)write(6,*)'READ_HDRAOB: no messages/reports ' + return + end if + write(6,*)'READ_HDRAOB: messages/reports = ',nmsg,'/',ntb,' ntread = ',ntread + + + if(qob .and. print_verbose) write(6,*)'READ_HDRAOB: time offset is ',toff,' hours.' +!------------------------------------------------------------------------ + + + allocate(list_stations(ntb)) +! loop over convinfo file entries; operate on matches + + allocate(cdata_all(nreal,maxobs)) + cdata_all=zero + nread=0 + ntest=0 + nvtest=0 + nchanl=0 + ilon=2 + ilat=3 + loop_convinfo: do nx=1,ntread + + call closbf(lunin) + open(lunin,file=infile,form='unformatted') + call openbf(lunin,'IN',lunin) + call datelen(10) + +! Big loop over hdraob file + + ntb = 0 + nmsg = 0 + disterrmax=-9999.0_r_kind + irec = 0 + loop_msg: do while (ireadmg(lunin,subset,idate)== 0) + irec = irec + 1 + if(irec < nrec_start) cycle loop_msg + + 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 + nc=tab(ntb,1) + kx = ictype(nc) + if(nc <= 0 .or. tab(ntb,2) /= kx) then + if(nc /= 0)write(6,*) obstype,nc,(tab(ntb,i),i=1,3),ntb + cycle loop_readsb + end if + +! Extract type, date, and location information + call ufbint(lunin,hdr,2,1,iret,hdstr) + igroup=hdr(1) + istation=hdr(2) + id=1000*igroup+istation +! write(6,*) id,igroup,istation,hdr(1),hdr(2) + if(igroup < 0 .or. istation < 0 .or. id >= 100000 .or. id <= 0)then + if(print_verbose)write(6,*) ' hdr ',hdr + cycle loop_readsb + end if + write(c_station_id,'(i5,3x)') id + call ufbint(lunin,levdat,2,maxlevs,levs,levstr) + if(levdat(1,2) > 1.e8 .and. (obstype == 'q' .or. obstype == 't'))cycle loop_readsb + call ufbint(lunin,obsdat,8,maxlevs,levs,obstr) + + if(psob)levs=1 + + if(qob)then + call ufbint(lunin,hdrtype,1,1,iret,hdtypestr) + r_station=hdrtype(1) + read(id_station,*) stntype + if(index(stntype,'IUD') /=0 .or. index( stntype,'IUX') /=0)then + if(kx == 119 .or. kx == 219 .or. kx == 218)then + write(6,*) ' inconsistent descent ',kx,id,levs,stntype + end if + else + if(kx == 118 .or. kx == 217)then + write(6,*) ' inconsistent ascent ',kx,id,levs,stntype + end if + end if + end if +! Combine height obs into single array and divide by g + do k=1,levs +! if(id == 47582)then +! write(6,*)id,levdat(1,k),k,(obsdat(i,k),i=1,8) +! end if + levdat(1,k)=.01*levdat(1,k) + if(levdat(2,k) > 1.e8)then + if(obsdat(4,k) < 1.e8)levdat(2,k)=obsdat(4,k)/grav + else + levdat(2,k)=levdat(2,k)/grav + end if + end do + +!------------------------------------------------------------------------ + + call ufbint(lunin,hdr2,8,1,iret,hdstr2) + dlat_earth_deg=hdr2(7) + dlon_earth_deg=hdr2(8) + if(abs(dlat_earth_deg)>r90 ) then + if(print_verbose)write(6,*) ' invalid lat ',id,dlat_earth_deg + cycle loop_readsb + end if + if(dlon_earth_deg >= r360)dlon_earth_deg=dlon_earth_deg-r360 + if(dlon_earth_deg < zero)dlon_earth_deg=dlon_earth_deg+r360 + dlon_earth=dlon_earth_deg*deg2rad + dlat_earth=dlat_earth_deg*deg2rad + 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 +! Interpolate guess 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 + pqm=0 + tqm=0 + qqm=0 + wqm=0 + if(levdat(1,2) < 1.e8)then + do k=1,levs + plevs(k)=0.1*levdat(1,k) ! convert mb to cb + end do + else +! If pressure is missing get an approximate pressure from guess field. This +! pressure is used for estimating ob error only. +! Note that balloon drift is not included in this calculation. + do kk=1,nsig + hgtl(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 + do k=1,levs + obheight=levdat(2,k) + if(obheight < hgtl(1))then + plevs(k)=presl(1) + else if (obheight > hgtl(nsig))then +! Reject ob by making pressure zero + plevs(k) = 0.0 + else + do kk=1,nsig-1 + if(hgtl(kk) <= obheight .and. hgtl(kk+1) > obheight)then + fact=(hgtl(kk+1)-obheight)/(hgtl(kk+1)-hgtl(kk)) + plevs(k)=exp(fact*log(presl(kk))+(1.-fact)*log(presl(kk+1))) + wqm(k)=1 + end if + end do + end if + end do + if(uvob .and. kx /= 218 )write(6,*) ' inconsistent kx ',kx + end if +! Compute depth of guess pressure layersat observation location + if (.not.twodvar_regional .and. levs > 1) then + do kk=1,nsig-1 + dpres(kk)=presl(kk)-presl(kk+1) + end do + endif + 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 minutes relative to historic date + rminan=minan + + 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 minutes 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 + + +! get observation launch time relative to analysis time + idate5(1)=hdr2(1) + idate5(2)=hdr2(2) + idate5(3)=hdr2(3) + idate5(4)=hdr2(4) + idate5(5)=hdr2(5) + call w3fs21(idate5,minobs) ! obs launch time in minutes relative to historic date + rminobs=minobs+hdr2(6)*r60inv + + timeobs=(rminobs-rminan)*r60inv + + t4dv=timeobs + toff + zeps=1.0e-8_r_kind + if (t4dv -zeps) t4dv=zero + if (t4dv>winlen.and.t4dv=10.0_r_kind) var_jb(1,k)=zero + enddo + endif + if (tob) then + ierr_t=0 + do i =1,maxsub_t + if( icsubtype(nc) == isuble_t(itypex,i) ) then + ierr_t=i+1 + exit + else if( i == maxsub_t .and. icsubtype(nc) /= isuble_t(itypex,i)) then + ncount_t=ncount_t+1 + do j=1,maxsub_t + if(isuble_t(itypex,j) ==0 ) then + ierr_t=j+1 + exit + endif + enddo + if( ncount_t ==1) then + write(6,*) 'READ_HDRAOB,WARNING!! tob:cannot find subtype in the error,& + table,itype=',itypex,icsubtype(nc) + write(6,*) 'read error table at column subtype as 0,error table column=',ierr_t + endif + endif + enddo + do k=1,levs + ppb=plevs(k)*10. + if(ppb <= zero .or. ppb > r2000)then + cycle loop_readsb + end if + if(ppb>=etabl_t(itypex,1,1)) k1_t=1 + do kl=1,32 + if(ppb>=etabl_t(itypex,kl+1,1).and.ppb<=etabl_t(itypex,kl,1)) k1_t=kl + end do + if(ppb<=etabl_t(itypex,33,1)) k1_t=5 + k2_t=k1_t+1 + ediff_t = etabl_t(itypex,k2_t,1)-etabl_t(itypex,k1_t,1) + if (abs(ediff_t) > tiny_r_kind) then + del_t = (ppb-etabl_t(itypex,k1_t,1))/ediff_t + else + del_t = huge_r_kind + endif + del_t=max(zero,min(del_t,one)) +! Temperature error + if(oberrflg)then + obserr(3,k)=(one-del_t)*etabl_t(itypex,k1_t,ierr_t)+del_t*etabl_t(itypex,k2_t,ierr_t) + obserr(3,k)=max(obserr(3,k),terrmin) + endif +!Temperature b + var_jb(3,k)=(one-del_t)*btabl_t(itypex,k1_t,ierr_t)+del_t*btabl_t(itypex,k2_t,ierr_t) + var_jb(3,k)=max(var_jb(3,k),zero) + if (var_jb(3,k) >=10.0_r_kind) var_jb(3,k)=zero + enddo + endif + if (qob) then + ierr_q=0 + do i =1,maxsub_q + if( icsubtype(nc) == isuble_q(itypex,i) ) then + ierr_q=i+1 + exit + else if( i == maxsub_q .and. icsubtype(nc) /= isuble_q(itypex,i)) then + ncount_q=ncount_q+1 + do j=1,maxsub_q + if(isuble_q(itypex,j) ==0 ) then + ierr_q=j+1 + exit + endif + enddo + if(ncount_q ==1 ) then + write(6,*) 'READ_HDRAOB,WARNING!! qob:cannot find subtype in the & + error table,itype=',itypex,icsubtype(nc) + write(6,*) 'read error table at column subtype as 0,error table column=',ierr_q + endif + endif + enddo + do k=1,levs + ppb=plevs(k)*10. + ppb=max(zero,min(ppb,r2000)) + if(ppb>=etabl_q(itypex,1,1)) k1_q=1 + do kl=1,32 + if(ppb>=etabl_q(itypex,kl+1,1).and.ppb<=etabl_q(itypex,kl,1)) k1_q=kl + end do + if(ppb<=etabl_q(itypex,33,1)) k1_q=5 + k2_q=k1_q+1 + ediff_q = etabl_q(itypex,k2_q,1)-etabl_q(itypex,k1_q,1) + if (abs(ediff_q) > tiny_r_kind) then + del_q = (ppb-etabl_q(itypex,k1_q,1))/ediff_q + else + del_q = huge_r_kind + endif + del_q=max(zero,min(del_q,one)) +! Humidity error + if(oberrflg)then + obserr(2,k)=(one-del_q)*etabl_q(itypex,k1_q,ierr_q)+del_q*etabl_q(itypex,k2_q,ierr_q) + obserr(2,k)=max(obserr(2,k),qerrmin) + endif +!Humidity b + var_jb(2,k)=(one-del_q)*btabl_q(itypex,k1_q,ierr_q)+del_q*btabl_q(itypex,k2_q,ierr_q) + var_jb(2,k)=max(var_jb(2,k),zero) + if (var_jb(2,k) >=10.0_r_kind) var_jb(2,k)=zero + enddo + endif + if (uvob) then + ierr_uv=0 + do i =1,maxsub_uv + if( icsubtype(nc) == isuble_uv(itypex,i) ) then + ierr_uv=i+1 + exit + else if( i == maxsub_uv .and. icsubtype(nc) /= isuble_uv(itypex,i)) then + ncount_uv=ncount_uv+1 + do j=1,maxsub_uv + if(isuble_uv(itypex,j) ==0 ) then + ierr_uv=j+1 + exit + endif + enddo + if( ncount_uv == 1) then + write(6,*) 'READ_HDRAOB,WARNING!! uvob:cannot find subtype in the error,& + table,itype=',itypex,icsubtype(nc) + write(6,*) 'read error table at column subtype as 0,error table column=',ierr_uv + endif + endif + enddo +! This has to be redone since no pressure for the ob. available. + do k=1,levs + ppb=plevs(k)*10. + ppb=max(zero,min(ppb,r2000)) + if(ppb>=etabl_uv(itypex,1,1)) k1_uv=1 + do kl=1,32 + if(ppb>=etabl_uv(itypex,kl+1,1).and.ppb<=etabl_uv(itypex,kl,1)) k1_uv=kl + end do + if(ppb<=etabl_uv(itypex,33,1)) k1_uv=5 + k2_uv=k1_uv+1 + ediff_uv = etabl_uv(itypex,k2_uv,1)-etabl_uv(itypex,k1_uv,1) + if (abs(ediff_uv) > tiny_r_kind) then + del_uv = (ppb-etabl_uv(itypex,k1_uv,1))/ediff_uv + else + del_uv = huge_r_kind + endif + del_uv=max(zero,min(del_uv,one)) +! Wind error + obserr(5,k)=(one-del_uv)*etabl_uv(itypex,k1_uv,ierr_uv)+del_uv*etabl_uv(itypex,k2_uv,ierr_uv) + obserr(5,k)=max(obserr(5,k),werrmin) +!Wind b + var_jb(5,k)=(one-del_uv)*btabl_uv(itypex,k1_uv,ierr_uv)+del_uv*btabl_uv(itypex,k2_uv,ierr_uv) + var_jb(5,k)=max(var_jb(5,k),zero) + if (var_jb(5,k) >=10.0_r_kind) var_jb(5,k)=zero + enddo + endif + else + do k=1,levs + itypex=kx + ppb=plevs(k)*10. + if(ppb <= zero .or. ppb > r2000)then + cycle loop_readsb + end if + if(ppb>=etabl(itypex,1,1)) k1=1 + do kl=1,32 + if(ppb>=etabl(itypex,kl+1,1).and.ppb<=etabl(itypex,kl,1)) k1=kl + end do + if(ppb<=etabl(itypex,33,1)) k1=5 + k2=k1+1 + if (abs(ediff) > tiny_r_kind) then + del = (ppb-etabl(itypex,k1,1))/ediff + else + del = huge_r_kind + endif + del=max(zero,min(del,one)) + obserr(3,k)=(one-del)*etabl(itypex,k1,2)+del*etabl(itypex,k2,2) + obserr(2,k)=(one-del)*etabl(itypex,k1,3)+del*etabl(itypex,k2,3) + obserr(5,k)=(one-del)*etabl(itypex,k1,4)+del*etabl(itypex,k2,4) + obserr(1,k)=(one-del)*etabl(itypex,k1,5)+del*etabl(itypex,k2,5) + obserr(7,k)=(one-del)*etabl(itypex,k1,6)+del*etabl(itypex,k2,6) + + obserr(3,k)=max(obserr(3,k),terrmin) + obserr(2,k)=max(obserr(2,k),qerrmin) + obserr(5,k)=max(obserr(5,k),werrmin) + obserr(1,k)=max(obserr(1,k),perrmin) + enddo + endif ! endif for njqc + + + +! If temperature ob, extract information regarding virtual +! versus sensible temperature + + + call ufbint(lunin,hdr3,1,1,iret,hdstr3) + stnelev=hdr3(1) + if(abs(hdr2(7)) > r90 .or. abs(hdr2(8)) > 720._r_kind)then + if(print_verbose)write(6,*) ' invalid lat,lon ',id,obstype,hdr2(7),hdr2(8) + cycle loop_readsb + end if + toocold=.false. + LOOP_K_LEVS: do k=1,levs + if(tob .or. qob .or. psob)then + if(levdat(1,k) < 1. .or. levdat(1,k) > 1500.)then + if(print_verbose)write(6,*) ' invalid pressure ',id,k,levs,levdat(1,k) + cycle LOOP_K_LEVS + end if + end if + if(obsdat(1,k) < 0. .or. obsdat(1,k) > 900000.) then + if(print_verbose)write(6,*) ' invalid change in time ',id,obstype,k,obsdat(1,k) + if(tob) tqm(k)=2 + if(qob) qqm(k)=2 + if(psob) pqm(k)=2 + if(uvob) wqm(k)=2 + obsdat(1,k) = 0. + end if + if(abs(obsdat(2,k)) < 10. .and. abs(obsdat(3,k)) < 10.) then + dlat_earth_deg=hdr2(7)+obsdat(2,k) + dlon_earth_deg=hdr2(8)+obsdat(3,k) + else + if(print_verbose)write(6,*) ' invalid change in lat/lon ',id,k,obstype,obsdat(2,k),obsdat(3,k) + if(tob) tqm(k)=2 + if(qob) qqm(k)=2 + if(psob) pqm(k)=2 + if(uvob) wqm(k)=2 + dlat_earth_deg=hdr2(7) + dlon_earth_deg=hdr2(8) + end if + if(dlon_earth_deg >= r360)dlon_earth_deg=dlon_earth_deg-r360 + if(dlon_earth_deg < zero)dlon_earth_deg=dlon_earth_deg+r360 + dlon_earth=dlon_earth_deg*deg2rad + dlat_earth=dlat_earth_deg*deg2rad + + 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 +! 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 + + + +! Set usage variable + usage = zero + if(icuse(nc) <= 0)usage=100._r_kind + if(pqm(k) >=lim_qm )usage=102._r_kind + + + if(plevs(k) < 0.0001_r_kind) then + write(*,*) 'warning: obs pressure is too small:',kx,k,plevs(k),id + cycle + endif + dlnpob=log(plevs(k)) ! ln(pressure in cb) + + call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz) + +! Extract pressure level and quality marks + +! Temperature + if(tob) then + if(obsdat(7,k) < 100. .or. obsdat(7,k) > 400.) then + if(print_verbose)write(6,*)id,'invalid temp',k,levs,obsdat(7,k),plevs(k) + tqm(k)=12 + usage=108._r_kind + cycle loop_k_levs + end if + ndata=ndata+1 + nodata=nodata+1 + iout=ndata + if(ndata > maxobs) then + write(6,*)'READ_HDRAOB: ***WARNING*** ndata > maxobs for ',obstype + ndata = maxobs + end if + ppb=plevs(k)*10. + if(levs > 100 .or. plevs(1)-plevs(levs) < .01)then + + call errormod_hdraob(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) + + else + + call errormod(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) + end if + toe=obserr(3,k)*errout + if(ppb < r100)toe=toe*r1_2 + cdata_all(1,iout)=toe ! temperature error + cdata_all(2,iout)=dlon ! grid relative longitude + cdata_all(3,iout)=dlat ! grid relative latitude + cdata_all(4,iout)=dlnpob ! ln(pressure in cb) + + cdata_all(5,iout)=obsdat(7,k) ! temperature ob. + cdata_all(6,iout)=rstation_id ! station id + cdata_all(7,iout)=t4dv+obsdat(1,k)*r60inv*r60inv ! time + cdata_all(8,iout)=nc ! type + qtflg=one + cdata_all(9,iout)=qtflg ! qtflg (virtual temperature flag) + cdata_all(10,iout)=tqm(k) ! quality mark + cdata_all(11,iout)=obserr(3,k) ! original obs error + cdata_all(12,iout)=usage ! usage parameter + cdata_all(13,iout)=idomsfc ! dominate surface type + cdata_all(14,iout)=tsavg ! skin temperature + cdata_all(15,iout)=ff10 ! 10 meter wind factor + cdata_all(16,iout)=sfcr ! surface roughness + cdata_all(17,iout)=dlon_earth_deg ! earth relative longitude (degrees) + cdata_all(18,iout)=dlat_earth_deg ! earth relative latitude (degrees) + cdata_all(19,iout)=stnelev ! station elevation (m) + cdata_all(20,iout)=levdat(2,k) ! observation height (m) + cdata_all(21,iout)=zz ! terrain height at ob location + cdata_all(22,iout)='88888888' ! provider name + cdata_all(23,iout)='HDRAOB' ! subprovider name + cdata_all(24,iout)=2 ! cat + cdata_all(25,iout)=var_jb(3,k) ! non linear qc for T + if(perturb_obs)cdata_all(nreal,iout)=ran01dom()*perturb_fact ! t perturbation + if (twodvar_regional) & + call adjust_error(cdata_all(17,iout),cdata_all(18,iout),cdata_all(11,iout),cdata_all(1,iout)) + if(kx == 119)then + if(usage < 100.)then + newstation=.true. + do i=1,nstations + if(list_stations(i) == id) then + newstation=.false. + exit + end if + end do + if(newstation)then + nstations=nstations+1 + list_stations(nstations)=id + end if + end if + end if + +! Winds + else if(uvob) then + if(obsdat(6,k) < 0. .or. obsdat(6,k) > 360.) then + wqm(k)=12 + usage=108._r_kind + if(print_verbose)write(6,*)id,'invalid dir ', k,levs,obsdat(5,k),obsdat(6,k),plevs(k) + cycle loop_k_levs + end if + if(obsdat(5,k) < 0. .or. obsdat(5,k) > 300.) then + wqm(k)=12 + usage=108._r_kind + if(print_verbose)write(6,*)id,'invalid spd ', k,levs,obsdat(5,k),obsdat(6,k),plevs(k) + cycle loop_k_levs + end if + ndata=ndata+1 + nodata=nodata+2 + iout=ndata + if(ndata > maxobs) then + write(6,*)'READ_HDRAOB: ***WARNING*** ndata > maxobs for ',obstype + ndata = maxobs + end if + if(levs > 100 .or. plevs(1)-plevs(levs) < .01)then + call errormod_hdraob(pqm,qqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) + else + + call errormod(pqm,wqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) + end if + woe=obserr(5,k)*errout + if(obsdat(1,k) < r50)woe=woe*r1_2 + oelev=levdat(2,k) + +! Rotate winds to rotated coordinate + uob = -obsdat(5,k)*sin(obsdat(6,k)*deg2rad) ! u-wind component + vob = -obsdat(5,k)*cos(obsdat(6,k)*deg2rad) ! v-wind component + + + if(regional .and. .not. fv3_regional)then + u0=uob + v0=vob + call rotate_wind_ll2xy(u0,v0,uob,vob,dlon_earth,dlon,dlat) + if(diagnostic_reg) then + call rotate_wind_xy2ll(uob,vob,u00,v00,dlon_earth,dlon,dlat) + nvtest=nvtest+1 + disterr=sqrt((u0-u00)**2+(v0-v00)**2) + vdisterrmax=max(vdisterrmax,disterr) + end if + endif + + cdata_all(1,iout)=woe ! wind error + cdata_all(2,iout)=dlon ! grid relative longitude + cdata_all(3,iout)=dlat ! grid relative latitude + cdata_all(4,iout)=dlnpob ! ln(pressure in cb) + cdata_all(5,iout)=oelev ! height of observation + cdata_all(6,iout)=uob ! u obs + cdata_all(7,iout)=vob ! v obs + cdata_all(8,iout)=rstation_id ! station id + cdata_all(9,iout)=t4dv+obsdat(1,k)*r60inv*r60inv ! time + cdata_all(10,iout)=nc ! type + cdata_all(11,iout)=stnelev ! station elevation + cdata_all(12,iout)=wqm(k) ! quality mark + cdata_all(13,iout)=obserr(5,k) ! original obs error + cdata_all(14,iout)=usage ! usage parameter + cdata_all(15,iout)=idomsfc ! dominate surface type + cdata_all(16,iout)=tsavg ! skin temperature + cdata_all(17,iout)=ff10 ! 10 meter wind factor + cdata_all(18,iout)=sfcr ! surface roughness + cdata_all(19,iout)=dlon_earth_deg ! earth relative longitude (degrees) + cdata_all(20,iout)=dlat_earth_deg ! earth relative latitude (degrees) + cdata_all(21,iout)=zz ! terrain height at ob location + cdata_all(22,iout)='88888888' ! provider name + cdata_all(23,iout)='HDRAOB' ! subprovider name + cdata_all(24,iout)=2 ! cat + cdata_all(25,iout)=var_jb(5,k) ! non linear qc parameter + cdata_all(26,iout)=one ! hilbert curve weight, modified later + if(perturb_obs)then + cdata_all(27,iout)=ran01dom()*perturb_fact ! u perturbation + cdata_all(28,iout)=ran01dom()*perturb_fact ! v perturbation + endif +! if(kx == 219 .and. k == 5)write(6,*) k,kx,c_station_id,(cdata_all(i,iout),i=1,25) + if(kx == 219 .or. kx == 218)then + if(usage < 100.)then + newstation=.true. + do i=1,nstations + if(list_stations(i) == id) then + newstation=.false. + exit + end if + end do + if(newstation)then + nstations=nstations+1 + list_stations(nstations)=id + end if + end if + end if + +! Specific humidity + else if(qob) then + if(obsdat(8,k) < 100. .or. obsdat(8,k) > 350.) then + if(print_verbose)write(6,*)id,'invalid td ', k,levs,obsdat(7,k),obsdat(8,k),plevs(k) + qqm(k)=12. + usage=108._r_kind + cycle loop_k_levs + end if + ndata=ndata+1 + nodata=nodata+1 + iout=ndata + if(ndata > maxobs) then + write(6,*)'READ_HDRAOB: ***WARNING*** ndata > maxobs for ',obstype + ndata = maxobs + end if + if(levs > 100 .or. plevs(1)-plevs(levs) < .01)then + call errormod_hdraob(pqm,qqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) + else + + call errormod(pqm,qqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) + end if + qoe=obserr(2,k)*errout*one_tenth + if(obsdat(7,k) < 215. .or. toocold) then + qqm(k)=12. + if(k > 1)toocold=.true. + usage=108._r_kind + end if +! Need to convert from td to q + call fpvsx_ad(obsdat(8,k),es,dummy,dummy,.false.) + qobcon = eps*es/(plevs(k)-omeps*es) + cdata_all(1,iout)=qoe ! q error + cdata_all(2,iout)=dlon ! grid relative longitude + cdata_all(3,iout)=dlat ! grid relative latitude + cdata_all(4,iout)=dlnpob ! ln(pressure in cb) + cdata_all(5,iout)=qobcon ! q ob + cdata_all(6,iout)=rstation_id ! station id + cdata_all(7,iout)=t4dv+obsdat(1,k)*r60inv*r60inv ! time + cdata_all(8,iout)=nc ! type + cdata_all(9,iout)=0.2_r_kind ! q max error + cdata_all(10,iout)=obsdat(7,k) ! dry temperature (obs is tv) + cdata_all(11,iout)=qqm(k) ! quality mark + cdata_all(12,iout)=obserr(2,k)*one_tenth ! original obs error + cdata_all(13,iout)=usage ! usage parameter + cdata_all(14,iout)=idomsfc ! dominate surface type + cdata_all(15,iout)=dlon_earth_deg ! earth relative longitude (degrees) + cdata_all(16,iout)=dlat_earth_deg ! earth relative latitude (degrees) + cdata_all(17,iout)=stnelev ! station elevation (m) + cdata_all(18,iout)=levdat(2,k) ! observation height (m) + cdata_all(19,iout)=zz ! terrain height at ob location + cdata_all(20,iout)='88888888' ! provider name + cdata_all(21,iout)='HDRAOB' ! subprovider name + cdata_all(22,iout)=2 ! cat + cdata_all(23,iout)=var_jb(2,k) ! non linear qc b parameter + if(perturb_obs)cdata_all(24,iout)=ran01dom()*perturb_fact ! q perturbation + if (twodvar_regional) & + call adjust_error(cdata_all(15,iout),cdata_all(16,iout),cdata_all(12,iout),cdata_all(1,iout)) + + if(kx == 119)then + if(usage < 100.)then + newstation=.true. + do i=1,nstations + if(list_stations(i) == id) then + newstation=.false. + exit + end if + end do + if(newstation)then + nstations=nstations+1 + list_stations(nstations)=id + end if + end if + end if + else if(psob) then + if(obsdat(7,k) > 400. .or. stnelev > 7000.) pqm(k)=12. + if(plevs(k) > 200. .or. plevs(k) < 50.) pqm(k)=12. + if(obsdat(4,k)/grav > 5000.)pqm(k)=12. + if(pqm(k) > lim_qm)then + usage=108._r_kind + exit loop_k_levs + end if + ndata=ndata+1 + nodata=nodata+1 + iout=ndata + if(ndata > maxobs) then + write(6,*)'READ_HDRAOB: ***WARNING*** ndata > maxobs for ',obstype + ndata = maxobs + end if + poe=obserr(1,k)*one_tenth ! convert from mb to cb + cdata_all(1,iout)=poe ! surface pressure error (cb) + cdata_all(2,iout)=dlon ! grid relative longitude + cdata_all(3,iout)=dlat ! grid relative latitude + + cdata_all(4,iout)=plevs(k) ! pressure (in cb) + + cdata_all(5,iout)=obsdat(4,k)/grav ! surface height + cdata_all(6,iout)=obsdat(7,k) ! surface temperature + cdata_all(7,iout)=rstation_id ! station id + cdata_all(8,iout)=t4dv ! time + cdata_all(9,iout)=nc ! type + cdata_all(10,iout)=pqm(k) ! quality mark + cdata_all(11,iout)=obserr(1,k)*one_tenth ! original obs error (cb) + cdata_all(12,iout)=usage ! usage parameter + cdata_all(13,iout)=idomsfc ! dominate surface type + 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)=stnelev ! station elevation (m) + cdata_all(17,iout)=zz ! terrain height at ob location + cdata_all(18,iout)='88888888' ! provider name + cdata_all(19,iout)='HDRAOB' ! subprovider name + cdata_all(20,iout)=var_jb(1,k) ! non linear qc b parameter + if(perturb_obs)cdata_all(21,iout)=ran01dom()*perturb_fact ! ps perturbation + if (twodvar_regional) & + call adjust_error(cdata_all(14,iout),cdata_all(15,iout),cdata_all(11,iout),cdata_all(1,iout)) + + if(usage < 100.)then + newstation=.true. + do i=1,nstations + if(list_stations(i) == id) then + newstation=.false. + exit + end if + end do + if(newstation)then + nstations=nstations+1 + list_stations(nstations)=id + end if + end if + end if + +! +! End k loop over levs + end do LOOP_K_LEVS + end do loop_readsb + +! +! End of bufr read loop + enddo loop_msg +! Close unit to bufr file + +! Normal exit + + enddo loop_convinfo! loops over convinfo entry matches + deallocate(lmsg,tab) + call closbf(lunin) + + if(print_verbose)write(6,*)'READ_HDRAOB: closbf(',lunin,')' + + close(lunin) + +! Write header record and data to output file for further processing + allocate(cdata_out(nreal,ndata)) + do i=1,ndata + do k=1,nreal + cdata_out(k,i)=cdata_all(k,i) + end do + end do + deallocate(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(diagnostic_reg .and. ntest>0) write(6,*)'READ_HDRAOB: ',& + 'ntest,disterrmax=',ntest,disterrmax + if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_HDRAOB: ',& + 'nvtest,vdisterrmax=',ntest,vdisterrmax + + if(tob)then + nhdt=nstations + if(nhdt > 0)allocate(hdtlist(nhdt)) + do i=1,nhdt + hdtlist(i)=list_stations(i) + end do + nodet=mype + write(6,*) ' number of high resolution t stations ',nstations + else if(qob)then + nhdq=nstations + if(nhdq > 0)allocate(hdqlist(nhdq)) + do i=1,nhdq + hdqlist(i)=list_stations(i) + end do + nodeq=mype + write(6,*) ' number of high resolution q stations ',nstations + else if(uvob)then + nhduv=nstations + if(nhduv > 0)allocate(hduvlist(nhduv)) + do i=1,nhduv + hduvlist(i)=list_stations(i) + end do + nodeuv=mype + write(6,*) ' number of high resolution uv stations ',nstations + else if(psob)then + nhdps=nstations + if(nhdps > 0)allocate(hdpslist(nhdps)) + do i=1,nhdps + hdpslist(i)=list_stations(i) + end do + nodeps=mype + write(6,*) ' number of high resolution ps stations ',nstations + end if + + deallocate(list_stations) +! End of routine + return + + end subroutine read_hdraob + +end module hdraobmod diff --git a/src/gsi/m_extOzone.F90 b/src/gsi/m_extOzone.F90 index 31cdedbebc..3d4b6783c1 100644 --- a/src/gsi/m_extOzone.F90 +++ b/src/gsi/m_extOzone.F90 @@ -1421,6 +1421,8 @@ subroutine ozlev_bufrread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, & call warn(myname_,' actual retained =',nodata) call warn(myname_,' size(ozout,2) =',maxobs) endif + call closbf(lunin) + close(lunin) ! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & ! ' obstype,nmrecs,ndata,nodata,no/ndata = ',dtype,nmrecs,ndata,nodata,real(nodata)/ndata diff --git a/src/gsi/mpimod.F90 b/src/gsi/mpimod.F90 index a249eadeac..7a4e597765 100644 --- a/src/gsi/mpimod.F90 +++ b/src/gsi/mpimod.F90 @@ -37,6 +37,7 @@ module mpimod use mpeu_mpif, only : mpi_real4 use mpeu_mpif, only : mpi_max use mpeu_mpif, only : mpi_min + use mpeu_mpif, only : mpi_character use mpeu_mpif, only : mpi_offset_kind use mpeu_mpif, only : mpi_info_null use mpeu_mpif, only : mpi_mode_rdonly @@ -109,6 +110,7 @@ module mpimod public :: mpi_mode_rdwr,mpi_byte public :: mpi_mode_create public :: mpi_mode_wronly + public :: mpi_character #ifdef HAVE_ESMF integer(i_kind) :: mpi_comm_world diff --git a/src/gsi/obs_para.f90 b/src/gsi/obs_para.f90 index e560fafe1f..530e946be6 100644 --- a/src/gsi/obs_para.f90 +++ b/src/gsi/obs_para.f90 @@ -104,7 +104,7 @@ subroutine obs_para(ndata,mype) if (dtype(is)=='lag') then ! lagrangian data call dislag(ndata(is,1),mm1,lunout,obsfile_all(is),dtype(is),& nobs_s) - nsat1(is)= nobs_sub(mm1,is) + nsat1(is)= nobs_sub(mm1,is) else obproc:do ii=1,npe if(nobs_sub(ii,is) > 0)then diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 294b6e1443..7fea9af263 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -84,6 +84,7 @@ module qcmod ! sub create_qcvars ! sub destroy_qcvars ! sub errormod +! sub errormod_hdraob ! sub errormod_aircraft ! sub setup_tzr_qc - set up QC with Tz retrieval ! sub tz_retrieval - Apply Tz retrieval @@ -162,6 +163,7 @@ module qcmod public :: create_qcvars public :: destroy_qcvars public :: errormod + public :: errormod_hdraob public :: errormod_aircraft public :: setup_tzr_qc public :: qc_ssmi @@ -634,10 +636,10 @@ subroutine errormod(pq,vq,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) implicit none integer(i_kind) ,intent(in ) :: levs,k,nsig,lim_qm - real(r_kind) ,dimension(255) ,intent(in ) :: plevs + real(r_kind) ,dimension(*) ,intent(in ) :: plevs real(r_kind) ,dimension(nsig) ,intent(in ) :: presl real(r_kind) ,dimension(nsig-1) ,intent(in ) :: dpres - integer(i_kind),dimension(255) ,intent(in ) :: pq,vq + integer(i_kind),dimension(*) ,intent(in ) :: pq,vq real(r_kind) ,intent(inout) :: errout integer(i_kind) n,l,ilev @@ -696,6 +698,117 @@ subroutine errormod(pq,vq,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) return end subroutine errormod + subroutine errormod_hdraob(pq,vq,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) +!$$$ subprogram documentation block +! . . . . +! subprogram: errormod +! prgmmr: derber org: np23 date: 2003-09-30 +! +! abstract: adjust observation error for conventional obs +! +! program history log: +! 2003-09-30 derber +! 2004-05-18 kleist, documentation +! 2004-10-26 kleist - add 0.5 half-layer factor +! 2006-02-15 treadon - add (l==levs,1) exit to upprof and dwprof loops +! 2006-12-20 Sienkiewicz multiply tiny_r_kind in errout div-by-zero +! check by expected largest value for numerator +! max(2*vmax) = max(dpres) ~= 5 cb +! 2008-04-23 safford - rm unused vars and uses +! 2008-09-05 lueken - merged ed's changes into q1fy09 code +! +! input argument list: +! pq - pressure quality mark +! vq - observation quality mark (t,q,wind) +! levs - number of levels in profile for observation +! plevs - observation pressures +! errout - observation error +! k - observation level +! presl - model pressure at half sigma levels +! dpres - delta pressure between model pressure levels +! nsig - number of vertical levels +! lim_qm - qc limit +! +! output argument list: +! errout - adjusted observation error +! +! attributes: +! language: f90 +! machine: ibm rs/6000 sp +! +!$$$ + implicit none + + integer(i_kind) ,intent(in ) :: levs,k,nsig,lim_qm + real(r_kind) ,dimension(*) ,intent(in ) :: plevs + real(r_kind) ,dimension(nsig) ,intent(in ) :: presl + real(r_kind) ,dimension(nsig-1) ,intent(in ) :: dpres + integer(i_kind),dimension(*) ,intent(in ) :: pq,vq + real(r_kind) ,intent(inout) :: errout + + integer(i_kind) n,l,ilev + real(r_kind):: vmag,pdiffu,pdiffd,con + + errout=one + if(levs == 1)return + ilev=1 + do n=2,nsig-1 + if(plevs(k) < presl(n))ilev=n + end do +! con=grav*500._r_kind/(273._r_kind*rd) +! vmag=min(max(half*dpres(ilev),r0_02*presl(ilev)),con*plevs(k)) + vmag=abs(dpres(ilev)) + +! vmag=max(half*dpres(ilev),r0_02*presl(ilev)) + pdiffu=vmag + pdiffd=vmag + if(pq(k) < lim_qm .and. vq(k) < lim_qm)then +! Move up through the profile. + l=k + +! Array plevs is only defined from l=1 to l=levs. Hence the check below + if (l+1<=levs) then + upprof: do while (abs(plevs(k)-plevs(l+1)) < vmag .and. l <= levs-1) + l=l+1 + if(pq(l) < lim_qm .and. vq(l) < lim_qm)then + pdiffu=abs(plevs(k)-plevs(l)) + exit upprof + end if + if (l==levs) exit upprof + end do upprof + endif + +! Reset the level and move down through the profile + l=k + +! The check (l>=2) ensures that plevs(l-1) is defined + if (l>=2) then + dwprof: do while (abs(plevs(l-1)-plevs(k)) < vmag .and. l >= 2) + l=l-1 + if(pq(l) < lim_qm .and. vq(l) < lim_qm)then + pdiffd=abs(plevs(l)-plevs(k)) + exit dwprof + end if + if (l==1) exit dwprof + end do dwprof + endif + +! Set adjusted error + if(plevs(k) > presl(1))then + errout=errout*(dpres(1)+plevs(k)-presl(1))/dpres(1) + else if(plevs(k) < presl(nsig))then + errout=errout*(dpres(nsig-1)+presl(nsig)-plevs(k))/dpres(nsig-1) + else + errout=sqrt(two*vmag/max(pdiffd+pdiffu,five*tiny_r_kind)) + end if + +! Quality marks indicate bad data. Set error to large value. + else + errout=1.e6_r_kind + end if + + return +end subroutine errormod_hdraob subroutine errormod_aircraft(pq,vq,levs,plevs,errout,k,presl,dpres,nsig,lim_qm,hdr3) !$$$ subprogram documentation block diff --git a/src/gsi/read_NASA_LaRC_cloud.f90 b/src/gsi/read_NASA_LaRC_cloud.f90 index cb6ceac65a..ec3f29e3fa 100644 --- a/src/gsi/read_NASA_LaRC_cloud.f90 +++ b/src/gsi/read_NASA_LaRC_cloud.f90 @@ -298,6 +298,7 @@ subroutine read_NASALaRC_cloud_bufr(satfile,atime,& enddo msg_report write(*,*) 'message/reports num=',nmsg,ntb call closbf(unit_in) + close(unit_in) numobs=ntb write(atime,'(I10)') idate @@ -409,6 +410,7 @@ subroutine read_NASALaRC_cloud_bufr_survey(satfile,east_time, west_time) enddo msg_report write(*,*) 'message/reports num=',nmsg,ntb call closbf(unit_in) + close(unit_in) write(*,'(2x,a10,a10,a11)') 'time_level','subset_num' DO i=1,num_obstime diff --git a/src/gsi/read_abi.f90 b/src/gsi/read_abi.f90 index 46c18a8516..5c4a5d1dbf 100644 --- a/src/gsi/read_abi.f90 +++ b/src/gsi/read_abi.f90 @@ -169,7 +169,7 @@ subroutine read_abi(mype,val_abi,ithin,rmesh,jsatid,& if (.not.assim) val_abi=zero ! Open bufr file. - call closbf(lnbufr) +! call closbf(lnbufr) open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) @@ -232,6 +232,7 @@ subroutine read_abi(mype,val_abi,ithin,rmesh,jsatid,& ! Reopen unit to bufr file call closbf(lnbufr) + close(lnbufr) open(lnbufr,file=infile,form='unformatted') call openbf(lnbufr,'IN',lnbufr) if(jsatid == 'gr' .or. jsatid == 'g16') kidsat = 270 @@ -494,6 +495,7 @@ subroutine read_abi(mype,val_abi,ithin,rmesh,jsatid,& enddo read_msg call closbf(lnbufr) + close(lnbufr) call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& nele,itxmax,nread,ndata,data_all,score_crit,nrec) diff --git a/src/gsi/read_ahi.f90 b/src/gsi/read_ahi.f90 index 3237f2e6fc..1d77a0ae2c 100644 --- a/src/gsi/read_ahi.f90 +++ b/src/gsi/read_ahi.f90 @@ -193,7 +193,7 @@ subroutine read_ahi(mype,val_img,ithin,rmesh,jsatid,gstime,& ! Open bufr file. - call closbf(lnbufr) +! call closbf(lnbufr) open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) @@ -219,6 +219,7 @@ subroutine read_ahi(mype,val_img,ithin,rmesh,jsatid,gstime,& allocate(data_all(nele,itxmax),nrec(itxmax)) call closbf(lnbufr) + close(lnbufr) open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) if(jsatid == 'himawari8') kidsat = 173 @@ -504,6 +505,7 @@ subroutine read_ahi(mype,val_img,ithin,rmesh,jsatid,gstime,& enddo read_loop enddo read_msg call closbf(lnbufr) + close(lnbufr) call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& nele,itxmax,nread,ndata,data_all,score_crit,nrec) diff --git a/src/gsi/read_airs.f90 b/src/gsi/read_airs.f90 index 62e4cc2ea6..ab6f770a91 100644 --- a/src/gsi/read_airs.f90 +++ b/src/gsi/read_airs.f90 @@ -855,6 +855,7 @@ subroutine read_airs(mype,val_airs,ithin,isfcalc,rmesh,jsatid,gstime,& enddo read_subset deallocate(allchan, chan_map, bufr_chan_test) call closbf(lnbufr) ! Close bufr file + close(lnbufr) ! If multiple tasks read input bufr file, allow each tasks to write out ! information it retained and then let single task merge files together diff --git a/src/gsi/read_amsr2.f90 b/src/gsi/read_amsr2.f90 index ed0429483b..a7b27abccc 100644 --- a/src/gsi/read_amsr2.f90 +++ b/src/gsi/read_amsr2.f90 @@ -428,6 +428,7 @@ subroutine read_amsr2(mype,val_amsr2,ithin,rmesh,jsatid,gstime,& enddo read_loop enddo call closbf(lnbufr) + close(lnbufr) num_obs=iobs-1 diff --git a/src/gsi/read_amsre.f90 b/src/gsi/read_amsre.f90 index 99a325c6c4..ef0c2ad2bb 100755 --- a/src/gsi/read_amsre.f90 +++ b/src/gsi/read_amsre.f90 @@ -645,6 +645,7 @@ subroutine read_amsre(mype,val_amsre,ithin,isfcalc,rmesh,jsatid,gstime,& enddo read_loop enddo read_msg call closbf(lnbufr) + close(lnbufr) ! If multiple tasks read input bufr file, allow each tasks to write out ! information it retained and then let single task merge files together diff --git a/src/gsi/read_anowbufr.f90 b/src/gsi/read_anowbufr.f90 index 6428e0d02e..e2b744eb6a 100644 --- a/src/gsi/read_anowbufr.f90 +++ b/src/gsi/read_anowbufr.f90 @@ -329,6 +329,7 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,& 1000 continue call closbf(lunin) + close(lunin) if(diagnostic_reg .and. & ntest > 0) write(6,*)'read_airnow_bufr: ',& diff --git a/src/gsi/read_atms.f90 b/src/gsi/read_atms.f90 index a0dce48f3a..1957e68109 100644 --- a/src/gsi/read_atms.f90 +++ b/src/gsi/read_atms.f90 @@ -374,7 +374,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,& end if ! Reopen unit to satellite bufr file - call closbf(lnbufr) +! call closbf(lnbufr) open(lnbufr,file=trim(infile2),form='unformatted',status = 'old', & iostat = ierr) if(ierr /= 0) cycle ears_db_loop @@ -492,6 +492,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,& end do read_loop end do read_subset call closbf(lnbufr) + close(lnbufr) end do ears_db_loop deallocate(data1b8) diff --git a/src/gsi/read_avhrr.f90 b/src/gsi/read_avhrr.f90 index ffc5a82052..8730b8bf13 100755 --- a/src/gsi/read_avhrr.f90 +++ b/src/gsi/read_avhrr.f90 @@ -532,6 +532,7 @@ subroutine read_avhrr(mype,val_avhrr,ithin,rmesh,jsatid,& enddo read_loop enddo read_msg call closbf(lnbufr) + close(lnbufr) call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& nele,itxmax,nread,ndata,data_all,score_crit,nrec) diff --git a/src/gsi/read_avhrr_navy.f90 b/src/gsi/read_avhrr_navy.f90 index 64a5b19907..dd5a64083a 100644 --- a/src/gsi/read_avhrr_navy.f90 +++ b/src/gsi/read_avhrr_navy.f90 @@ -496,6 +496,7 @@ subroutine read_avhrr_navy(mype,val_avhrr,ithin,rmesh,jsatid,& 900 continue call destroygrids call closbf(lnbufr) + close(lnbufr) if(diagnostic_reg.and.ntest>0) write(6,*)'READ_AVHRR_NAVY: ',& 'mype,ntest,disterrmax=',mype,ntest,disterrmax diff --git a/src/gsi/read_bufrtovs.f90 b/src/gsi/read_bufrtovs.f90 index 9e9795b009..cdba37f9e0 100644 --- a/src/gsi/read_bufrtovs.f90 +++ b/src/gsi/read_bufrtovs.f90 @@ -501,7 +501,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& end if ! Reopen unit to satellite bufr file - call closbf(lnbufr) +! call closbf(lnbufr) open(lnbufr,file=trim(infile2),form='unformatted',status = 'old',iostat=ierr) if(ierr /= 0) cycle ears_db_loop @@ -946,6 +946,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& enddo read_loop enddo read_subset call closbf(lnbufr) + close(lnbufr) if(llll > 1 .and. (amsua .or. amsub .or. mhs))then deallocate(data1b8x) diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index 91bac38c65..29067b3f69 100644 --- a/src/gsi/read_cris.f90 +++ b/src/gsi/read_cris.f90 @@ -404,7 +404,6 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& end if ! Open BUFR file - call closbf(lnbufr) open(lnbufr,file=trim(infile2),form='unformatted',status='old',iostat=ierr) if(ierr /= 0) cycle ears_db_loop @@ -843,6 +842,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& enddo read_subset call closbf(lnbufr) + close(lnbufr) end do ears_db_loop diff --git a/src/gsi/read_fl_hdob.f90 b/src/gsi/read_fl_hdob.f90 index bd2e9929cb..b3f4b3ba65 100644 --- a/src/gsi/read_fl_hdob.f90 +++ b/src/gsi/read_fl_hdob.f90 @@ -358,7 +358,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si ! Go through the bufr file to find out how mant subsets to process nmsg = 0 maxobs = 0 - call closbf(lunin) +! call closbf(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -371,6 +371,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si end do loop_readsb1 end do loop_msg1 call closbf(lunin) + close(lunin) write(6,*) 'READ_FL_HDOB: total number of data found in the bufr file ',maxobs,obstype write(6,*) 'READ_FL_HDOB: time offset is ',toff,' hours' @@ -391,7 +392,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si ilat = 3 ! Open bufr file again for reading - call closbf(lunin) +! call closbf(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -1158,6 +1159,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si ! Close unit to bufr file call closbf(lunin) + close(lunin) ! Deallocate arrays used for thinning data if (.not.use_all) then deallocate(presl_thin) @@ -1185,13 +1187,11 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si 'nvtest,vdisterrmax=',ntest,vdisterrmax if (ndata == 0) then - call closbf(lunin) write(6,*)'READ_FL_HDOB: no data to process' endif write(6,*)'READ_FL_HDOB: nreal=',nreal write(6,*)'READ_FL_HDOB: ntb,nread,ndata,nodata=',ntb,nread,ndata,nodata - close(lunin) ! End of routine return diff --git a/src/gsi/read_gmi.f90 b/src/gsi/read_gmi.f90 index 75b56b0e14..4361e7144a 100644 --- a/src/gsi/read_gmi.f90 +++ b/src/gsi/read_gmi.f90 @@ -501,6 +501,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& end do read_subset 690 continue call closbf(lnbufr) + close(lnbufr) num_obs=iobs-1 if( mype_sub==mype_root) write(6,*) 'READ_GMI: do_noise_reduction=', do_noise_reduction diff --git a/src/gsi/read_goesglm.f90 b/src/gsi/read_goesglm.f90 index f72e53623d..b2b77fb693 100644 --- a/src/gsi/read_goesglm.f90 +++ b/src/gsi/read_goesglm.f90 @@ -128,7 +128,7 @@ subroutine read_goesglm(nread,ndata,nodata,infile,obstype,lunout,twindin,sis) ! Open, then read date from BUFR file - call closbf(lunin) +! call closbf(lunin) open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -401,15 +401,12 @@ subroutine read_goesglm(nread,ndata,nodata,infile,obstype,lunout,twindin,sis) 'nvtest,vdisterrmax=',ntest,vdisterrmax if (ndata == 0) then - call closbf(lunin) - write(6,*)'READ_GOESGLM: closbf(',lunin,')' + write(6,*)'READ_GOESGLM: closbf(',lunin,') no data' endif close(lunin) - close(55) - ! End of routine return diff --git a/src/gsi/read_goesimg.f90 b/src/gsi/read_goesimg.f90 index e3ceaa7929..bf40a1f163 100644 --- a/src/gsi/read_goesimg.f90 +++ b/src/gsi/read_goesimg.f90 @@ -199,7 +199,6 @@ subroutine read_goesimg(mype,val_img,ithin,rmesh,jsatid,gstime,& ! Open bufr file. - call closbf(lnbufr) open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) @@ -444,6 +443,7 @@ subroutine read_goesimg(mype,val_img,ithin,rmesh,jsatid,gstime,& 900 continue call destroygrids call closbf(lnbufr) + close(lnbufr) if(diagnostic_reg.and.ntest>0) write(6,*)'READ_GOESIMG: ',& 'mype,ntest,disterrmax=',mype,ntest,disterrmax diff --git a/src/gsi/read_goesimgr_skycover.f90 b/src/gsi/read_goesimgr_skycover.f90 index e3c5432d82..dda9aad6f4 100644 --- a/src/gsi/read_goesimgr_skycover.f90 +++ b/src/gsi/read_goesimgr_skycover.f90 @@ -205,11 +205,11 @@ subroutine read_goesimgr_skycover(nread,ndata,nodata,infile,obstype,lunout,gsti ilat=3 ntb=0 - close(lunin) - call closbf(lunin) - open(lunin,file=trim(infile),form='unformatted') - call openbf(lunin,'IN',lunin) - call datelen(10) + 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) @@ -418,7 +418,7 @@ subroutine read_goesimgr_skycover(nread,ndata,nodata,infile,obstype,lunout,gsti deallocate(cdata_out) if (ndata == 0) then - write(6,*)myname,': closbf(',lunin,')' + write(6,*)myname,'no read_goesimgr_skycover data' endif close(lunin) diff --git a/src/gsi/read_goesndr.f90 b/src/gsi/read_goesndr.f90 index 8335b82d93..86fc1f0a5c 100644 --- a/src/gsi/read_goesndr.f90 +++ b/src/gsi/read_goesndr.f90 @@ -513,6 +513,7 @@ subroutine read_goesndr(mype,val_goes,ithin,rmesh,jsatid,infile,& end do read_loop end do read_subset call closbf(lnbufr) + close(lnbufr) ! If multiple tasks read input bufr file, allow each tasks to write out diff --git a/src/gsi/read_gps.f90 b/src/gsi/read_gps.f90 index fc3888dcf1..3d8379ee3b 100644 --- a/src/gsi/read_gps.f90 +++ b/src/gsi/read_gps.f90 @@ -460,6 +460,7 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & ! Close unit to input file call closbf(lnbufr) + close(lnbufr) nprof_gps = nmrecs write(6,*)'READ_GPS: # bad or missing data=', notgood diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90 index 0079a2c239..dc84e4bd64 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -414,7 +414,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& end if ! Open BUFR file - call closbf(lnbufr) +! call closbf(lnbufr) open(lnbufr,file=trim(infile2),form='unformatted',status='old',iostat=ierr) if(ierr /= 0) cycle ears_db_loop @@ -808,6 +808,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& enddo read_subset call closbf(lnbufr) + close(lnbufr) end do ears_db_loop diff --git a/src/gsi/read_l2bufr_mod.f90 b/src/gsi/read_l2bufr_mod.f90 index a124b75d3d..cdcb389ecf 100644 --- a/src/gsi/read_l2bufr_mod.f90 +++ b/src/gsi/read_l2bufr_mod.f90 @@ -395,6 +395,7 @@ subroutine radar_bufr_read_all(npe,mype) if(rite) write(6,*)'RADAR_BUFR_READ_ALL: NO RADARS KEPT IN radar_bufr_read_all, ',& 'continue without level 2 data' call closbf(inbufr) + close(inbufr) return end if call mpi_reduce(num_radars,num_radars_min,1,mpi_integer4,mpi_min,0,mpi_comm_world,ierror) @@ -496,6 +497,7 @@ subroutine radar_bufr_read_all(npe,mype) end if ! reopen and reread the file for data this time call closbf(inbufr) + close(inbufr) open(inbufr,file=infile,form='unformatted') call openbf(inbufr,'IN',inbufr) allocate(bins(6,nthisrad,num_radars_0),ibins(nthisrad,num_radars_0)) @@ -684,6 +686,7 @@ subroutine radar_bufr_read_all(npe,mype) end do ! end do while end do ! loop over blocks call closbf(inbufr) + close(inbufr) if (.not. allocated(ibins2)) allocate(ibins2(nthisrad,num_radars_0)) ibins2=0 @@ -931,6 +934,7 @@ subroutine radar_bufr_read_all(npe,mype) write(6,*)' deldistmin,maxall=',deldistminall,deldistmaxall end if close(inbufr) + close(inbufr) end if deallocate(bins_work,bins,ibins2) if(l2superob_only) then diff --git a/src/gsi/read_lidar.f90 b/src/gsi/read_lidar.f90 index 3112606eed..6d74de0802 100644 --- a/src/gsi/read_lidar.f90 +++ b/src/gsi/read_lidar.f90 @@ -305,6 +305,7 @@ subroutine read_lidar(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) ! Close unit to bufr file deallocate(cdata_all) call closbf(lunin) + close(lunin) ! End of routine return diff --git a/src/gsi/read_modsbufr.f90 b/src/gsi/read_modsbufr.f90 index db551844ff..f2fab23a86 100644 --- a/src/gsi/read_modsbufr.f90 +++ b/src/gsi/read_modsbufr.f90 @@ -573,6 +573,7 @@ subroutine read_modsbufr(nread,ndata,nodata,gstime,infile,obstype,lunout, & 1020 continue if (oberrflg) deallocate(etabl) call closbf(lunin) + close(lunin) if(regional)then if(diagnostic_reg.and.ntest > 0) write(6,*)'READ_MODSBUFR: ',& diff --git a/src/gsi/read_nasa_larc.f90 b/src/gsi/read_nasa_larc.f90 index 71c6e939a4..8dbbbb2e68 100644 --- a/src/gsi/read_nasa_larc.f90 +++ b/src/gsi/read_nasa_larc.f90 @@ -199,6 +199,7 @@ subroutine read_nasa_larc(nread,ndata,infile,obstype,lunout,twind,sis,nobs) endif ! call closbf(lunin) + close(lunin) return 200 continue write(6,*) 'read_nasa_larc, Warning : cannot find LaRC data file' diff --git a/src/gsi/read_nsstbufr.f90 b/src/gsi/read_nsstbufr.f90 index 31435280dc..2b8dc239ae 100644 --- a/src/gsi/read_nsstbufr.f90 +++ b/src/gsi/read_nsstbufr.f90 @@ -604,6 +604,7 @@ subroutine read_nsstbufr(nread,ndata,nodata,gstime,infile,obstype,lunout, & 1020 continue if (oberrflg) deallocate(etabl) call closbf(lunin) + close(lunin) if(diagnostic_reg.and.ntest > 0) write(6,*)'READ_NSSTBUFR: ',& 'ntest,disterrmax=',ntest,disterrmax diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 2c485a3ee6..076149c3fe 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -60,6 +60,7 @@ subroutine gsi_inquire (lbytes,lexist,filename,mype) ! machine: Linux-cluster ! !$$$ end documentation block +! use kinds, only: i_kind,i_llong @@ -326,6 +327,7 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) end if call closbf(lnbufr) + close(lnbufr) open(lnbufr,file=trim(filename),form='unformatted',status ='unknown') call openbf(lnbufr,'IN',lnbufr) call datelen(10) @@ -551,6 +553,7 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) end if call closbf(lnbufr) + close(lnbufr) end if if(lexist)then write(6,*)'read_obs_check: bufr file date is ',idate,trim(filename),' ',dtype,jsatid @@ -704,12 +707,13 @@ subroutine read_obs(ndata,mype) reduce_diag,nobs_sub,dval_use,hurricane_radar,l2rwthin use gsi_nstcouplermod, only: nst_gsi ! use gsi_nstcouplermod, only: gsi_nstcoupler_set + use hdraobmod, only: read_hdraob,nhdt,nhdq,nhduv,nhdps,hdtlist,hdqlist,hduvlist,hdpslist,nodet,nodeq,nodeuv,nodeps use qcmod, only: njqc,vadwnd_l2rw_qc,nvqc use gsi_4dvar, only: l4dvar use satthin, only: super_val,super_val1,superp,makegvals,getsfc,destroy_sfc - use mpimod, only: ierror,mpi_comm_world,mpi_sum,mpi_rtype,mpi_integer,npe,& + use mpimod, only: ierror,mpi_comm_world,mpi_sum,mpi_max,mpi_rtype,mpi_integer,npe,& setcomm - use constants, only: one,zero + use constants, only: one,zero,izero use converr, only: converr_read use converr_ps, only: converr_ps_read use converr_q, only: converr_q_read @@ -767,6 +771,7 @@ subroutine read_obs(ndata,mype) integer(i_kind):: npetot,npeextra,mmdat,nodata integer(i_kind):: iworld,iworld_group,next_mype,mm1,iix integer(i_kind):: mype_root + integer(i_kind),dimension(8):: nhd,nhd1 integer(i_kind):: minuse,lunsave,maxproc,minproc integer(i_kind),dimension(ndat):: npe_sub,npe_sub3,mpi_comm_sub,mype_root_sub,npe_order integer(i_kind),dimension(ndat):: ntasks1,ntasks @@ -836,7 +841,15 @@ subroutine read_obs(ndata,mype) deallocate(nrnd) endif - +! Set number of high definition stations to zero + nhdt=izero + nhdq=izero + nhduv=izero + nhdps=izero + nodet=izero + nodeq=izero + nodeuv=izero + nodeps=izero ! Set data class and number of reader tasks. Set logical flag to indicate ! type type of GPS data (if present) @@ -1252,10 +1265,15 @@ subroutine read_obs(ndata,mype) end if end do end if - if(obstype == 'rw')then + if(dfile(i) == 'uprair')then + use_prsl_full=.true. + use_hgtl_full = .true. + if(belong(i))use_hgtl_full_proc=.true. + if(belong(i))use_prsl_full_proc=.true. + else if(obstype == 'rw')then use_hgtl_full=.true. if(belong(i))use_hgtl_full_proc=.true. - else if(obstype == 'dbz')then + else if(obstype == 'dbz')then use_hgtl_full=.true. if(belong(i))use_hgtl_full_proc=.true. end if @@ -1389,6 +1407,10 @@ subroutine read_obs(ndata,mype) call read_fl_hdob(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,& prsl_full,nobs_sub1(1,i)) string='READ_FL_HDOB' + else if (index(infile,'uprair') /=0)then + call read_hdraob(nread,npuse,nouse,infile,obstype,lunout,twind,sis,& + prsl_full,hgtl_full,nobs_sub1(1,i),read_rec(i)) + string='READ_UPRAIR' else call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,& prsl_full,nobs_sub1(1,i),read_rec(i)) @@ -1452,6 +1474,11 @@ subroutine read_obs(ndata,mype) call read_satwnd(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,& prsl_full,nobs_sub1(1,i)) string='READ_SATWND' +! Process high resolution radiosonde data + else if (index(infile,'uprair') /=0)then + call read_hdraob(nread,npuse,nouse,infile,obstype,lunout,twind,sis,& + prsl_full,hgtl_full,nobs_sub1(1,i),read_rec(i)) + string='READ_UPRAIR' ! Process oscat winds which seperate from prepbufr elseif ( index(infile,'oscatbufr') /=0 ) then call read_sfcwnd(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,& @@ -1466,6 +1493,10 @@ subroutine read_obs(ndata,mype) call read_fl_hdob(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,& prsl_full,nobs_sub1(1,i)) string='READ_FL_HDOB' + else if (index(infile,'uprair') /=0)then + call read_hdraob(nread,npuse,nouse,infile,obstype,lunout,twind,sis,& + prsl_full,hgtl_full,nobs_sub1(1,i),read_rec(i)) + string='READ_UPRAIR' else call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,& prsl_full,nobs_sub1(1,i),read_rec(i)) @@ -1837,8 +1868,6 @@ subroutine read_obs(ndata,mype) end if ditype_select -! Close unit to data file - ! Accumulate data counts on "root" task if (mype_sub(mm1,i)==mype_root) then close(lunout) @@ -1891,6 +1920,43 @@ subroutine read_obs(ndata,mype) end if super_val1(0)=one deallocate(super_val) + nhd(1)=nhdt + nhd(2)=nhdq + nhd(3)=nhduv + nhd(4)=nhdps + nhd(5)=nodet + nhd(6)=nodeq + nhd(7)=nodeuv + nhd(8)=nodeps +! get number of high resolution stations on every processor + call mpi_allreduce(nhd,nhd1,8,mpi_integer,mpi_max,mpi_comm_world,ierror) + nhdt=nhd1(1) + nhdq=nhd1(2) + nhduv=nhd1(3) + nhdps=nhd1(4) + nodet=nhd1(5) + nodeq=nhd1(6) + nodeuv=nhd1(7) + nodeps=nhd1(8) + if(nhdt > 0)then + if(.not. allocated(hdtlist))allocate(hdtlist(nhdt)) + call mpi_bcast(hdtlist,nhdt,mpi_integer,nodet,mpi_comm_world,ierror) + end if + if(nhdq > 0) then + if(.not. allocated(hdqlist))allocate(hdqlist(nhdq)) + call mpi_bcast(hdqlist,nhdq,mpi_integer,nodeq,mpi_comm_world,ierror) + end if + if(nhduv > 0) then + if(.not. allocated(hduvlist))allocate(hduvlist(nhduv)) + call mpi_bcast(hduvlist,nhduv,mpi_integer,nodeuv,mpi_comm_world,ierror) + end if + if(nhdps > 0)then + if(.not. allocated(hdpslist))allocate(hdpslist(nhduv)) + call mpi_bcast(hdpslist,nhdps,mpi_integer,nodeps,mpi_comm_world,ierror) + end if + + if(mype == 0)write(6,*)'number of HD stations',nhdt,nhdq,nhduv,nhdps + if(mype == 0)write(6,*)'processors ',nodet,nodeq,nodeuv,nodeps ! Collect number of gps profiles (needed later for qc) call mpi_allreduce(nprof_gps1,nprof_gps,1,mpi_integer,mpi_sum,mpi_comm_world,ierror) diff --git a/src/gsi/read_ozone.f90 b/src/gsi/read_ozone.f90 index 74c95e09db..9e1cd94168 100644 --- a/src/gsi/read_ozone.f90 +++ b/src/gsi/read_ozone.f90 @@ -424,6 +424,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ! Loop back to read next profile end do read_loop1 + call closbf(lunin) + close(lunin) ! End of bufr ozone block ! Process GOME-2 data @@ -583,6 +585,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & end do obsloop + call closbf(lunin) + close(lunin) ! End of GOME bufr block @@ -753,6 +757,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ! End of loop over observations end do read_loop2 + call closbf(lunin) + close(lunin) ! End of OMI/OMPS-NM(or TC8) block @@ -862,6 +868,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ! Reopen unit to bufr file call closbf(lunin) + close(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -1056,6 +1063,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & end do read_loop4 + call closbf(lunin) + close(lunin) ! End of MLS bufr loop !Process OMPS LP data @@ -1222,6 +1231,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & enddo enddo read_loop5 + call closbf(lunin) + close(lunin) ! end of OMPS LP bufr loop endif @@ -1270,10 +1281,6 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & if(allocated(usage1))deallocate(usage1) endif -! Close unit to input data file - call closbf(lunin) - close(lunin) - ! Deallocate satthin arrays if (obstype == 'omi' .or. obstype == 'gome' .or. obstype=='ompsnm' .or. obstype == 'ompstc8' )call destroygrids diff --git a/src/gsi/read_pblh.f90 b/src/gsi/read_pblh.f90 index ff7d2f220b..a7b7c066a2 100644 --- a/src/gsi/read_pblh.f90 +++ b/src/gsi/read_pblh.f90 @@ -99,7 +99,7 @@ subroutine read_pblh(nread,ndata,nodata,infile,obstype,lunout,twindin,& ntest=0 nrtmax=0 ! # rpts to print per msg type (0=all) - call closbf(lunin) +! call closbf(lunin) open(lunin,file=trim(infile),form='unformatted') call mesgbc(lunin,msgt,icomp) call openbf(lunin,'IN',lunin) @@ -157,6 +157,7 @@ subroutine read_pblh(nread,ndata,nodata,infile,obstype,lunout,twindin,& ilon=2 ilat=3 call closbf(lunin) + close(lunin) open(lunin,file=trim(infile),form='unformatted') call mesgbc(lunin,msgt,icomp) call openbf(lunin,'IN',lunin) @@ -477,6 +478,7 @@ subroutine read_pblh(nread,ndata,nodata,infile,obstype,lunout,twindin,& write(*,*) ! closing linefeed, debug? call closbf(lunin) + close(lunin) ! Normal exit ! Write observation to scratch file @@ -486,9 +488,7 @@ subroutine read_pblh(nread,ndata,nodata,infile,obstype,lunout,twindin,& deallocate(cdata_all) if (ndata == 0) then - call closbf(lunin) - write(6,*)'READ_PREPFITS: closbf(',lunin,')' + write(6,*)'READ_PREPFITS no data' endif - close(lunin) end subroutine read_pblh diff --git a/src/gsi/read_pcp.f90 b/src/gsi/read_pcp.f90 index 929a608ee7..fe5fe456e6 100644 --- a/src/gsi/read_pcp.f90 +++ b/src/gsi/read_pcp.f90 @@ -147,7 +147,7 @@ subroutine read_pcp(nread,ndata,nodata,gstime,infile,lunout,obstype, & ! Open and read the bufr data - call closbf(lnbufr) +! call closbf(lnbufr) open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) @@ -354,6 +354,7 @@ subroutine read_pcp(nread,ndata,nodata,gstime,infile,lunout,obstype, & ! Jump here if there is a problem opening the bufr file 110 continue call closbf(lnbufr) + close(lnbufr) ! End of routine return diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index 9041d9702b..b22a7bf443 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -615,7 +615,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ncount_ps=0;ncount_q=0;ncount_t=0;ncount_uv=0;ncount_pw=0 ! Open, then read date from bufr data - call closbf(lunin) +! call closbf(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -910,6 +910,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& call closbf(lunin) + close(lunin) open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -1108,6 +1109,12 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if (sfctype) then call ufbint(lunin,r_prvstg,1,1,iret,prvstr) call ufbint(lunin,r_sprvstg,1,1,iret,sprvstr) + else if(kx == 120 .and. tob .or. qob .or. psob)then + c_prvstg=cspval + c_sprvstg='PREP' + else if(kx == 220 .and. uvob)then + c_prvstg=cspval + c_sprvstg='PREP' else c_prvstg=cspval c_sprvstg=cspval @@ -2343,7 +2350,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if(perturb_obs)cdata_all(24,iout)=ran01dom()*perturb_fact ! q perturbation if (twodvar_regional) & call adjust_error(cdata_all(15,iout),cdata_all(16,iout),cdata_all(12,iout),cdata_all(1,iout)) - + ! Total precipitable water (ssm/i) else if(pwob) then @@ -2908,8 +2915,6 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! ! End of bufr read loop enddo loop_msg -! Close unit to bufr file - call closbf(lunin) ! Deallocate arrays used for thinning data if (.not.use_all) then @@ -2926,6 +2931,9 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& enddo loop_convinfo! loops over convinfo entry matches deallocate(lmsg,tab,nrep) +! Close unit to bufr file + call closbf(lunin) + close(lunin) ! Apply hilbert curve for cross validation if requested @@ -3152,10 +3160,8 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_PREPBUFR: ',& 'nvtest,vdisterrmax=',ntest,vdisterrmax - call closbf(lunin) if(print_verbose)write(6,*)'READ_PREPBUFR: closbf(',lunin,')' - close(lunin) ! End of routine return diff --git a/src/gsi/read_radar.f90 b/src/gsi/read_radar.f90 index 7493e5d71c..2d7c039995 100644 --- a/src/gsi/read_radar.f90 +++ b/src/gsi/read_radar.f90 @@ -344,6 +344,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu end do loop0 call closbf(lnbufr) + close(lnbufr) ! enddo msg_report @@ -471,6 +472,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu ! Normal exit end if call closbf(lnbufr) + close(lnbufr) ! Print vadwnd table @@ -813,7 +815,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu end do - close(lnbufr) ! A simple unformatted fortran file should not be mixed with a bufr I/O + close(lnbufr) ! A simple unformatted fortran file should not be mixed with a bufr I/O LEVEL_TWO_READ_2: if(loop==0 .and. sis=='l2rw') then write(6,*)'READ_RADAR: ',trim(outmessage),' reached eof on 2/2.5/3 superob radar file' @@ -1249,6 +1251,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu ! Close unit to bufr file call closbf(lnbufr) + close(lnbufr) write(6,*)'READ_RADAR: ',trim(outmessage),' reached eof on 2.5/3 superob radar file.' @@ -1688,6 +1691,7 @@ subroutine read_radar(nread,ndata,nodata,infile,lunout,obstype,twind,sis,hgtl_fu close(25) end do ! end of loop, reading TDR so data files + close(lnbufr) else diff --git a/src/gsi/read_radarref_mosaic.f90 b/src/gsi/read_radarref_mosaic.f90 index 55e141d8b3..dab4da6e83 100644 --- a/src/gsi/read_radarref_mosaic.f90 +++ b/src/gsi/read_radarref_mosaic.f90 @@ -212,6 +212,7 @@ subroutine read_radarref_mosaic(nread,ndata,infile,obstype,lunout,twind,sis,nobs endif call closbf(lunin) + close(lunin) return 200 continue write(6,*) 'read_radarref_mosaic, Warning : cannot find radar data file' diff --git a/src/gsi/read_rapidscat.f90 b/src/gsi/read_rapidscat.f90 index dc48991940..c952383df0 100644 --- a/src/gsi/read_rapidscat.f90 +++ b/src/gsi/read_rapidscat.f90 @@ -244,7 +244,6 @@ subroutine read_rapidscat(nread,ndata,nodata,infile,obstype,lunout,gstime,twind, !! go through the satedump to find out how many subset to process !** Open and read data from bufr data file - call closbf(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -389,6 +388,7 @@ subroutine read_rapidscat(nread,ndata,nodata,infile,obstype,lunout,gstime,twind, endif call closbf(lunin) + close(lunin) open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -663,8 +663,6 @@ subroutine read_rapidscat(nread,ndata,nodata,infile,obstype,lunout,gstime,twind, 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) @@ -674,6 +672,7 @@ subroutine read_rapidscat(nread,ndata,nodata,infile,obstype,lunout,gstime,twind, enddo loop_convinfo! loops over convinfo entry matches deallocate(lmsg,tab,nrep) + call closbf(lunin) ! Write header record and data to output file for further processing allocate(iloc(ndata)) @@ -710,10 +709,7 @@ subroutine read_rapidscat(nread,ndata,nodata,infile,obstype,lunout,gstime,twind, if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_RAPIDSCAT: ',& 'nvtest,vdisterrmax=',ntest,vdisterrmax - if (ndata == 0) then - call closbf(lunin) - write(6,*)'READ_RAPIDSCAT: closbf(',lunin,')' - endif + write(6,*)'READ_RAPIDSCAT: closbf(',lunin,')' write(6,*) 'READ_RAPIDSCAT,nread,ndata,nreal,nodata=',nread,ndata,nreal,nodata diff --git a/src/gsi/read_saphir.f90 b/src/gsi/read_saphir.f90 index e4d7c55bde..3547f2cf5c 100644 --- a/src/gsi/read_saphir.f90 +++ b/src/gsi/read_saphir.f90 @@ -283,7 +283,6 @@ subroutine read_saphir(mype,val_tovs,ithin,isfcalc,& ! Reopen unit to satellite bufr file iob=1 - call closbf(lnbufr) open(lnbufr,file=trim(infile),form='unformatted',status = 'old',err = 500) call openbf(lnbufr,'IN',lnbufr) @@ -388,6 +387,7 @@ subroutine read_saphir(mype,val_tovs,ithin,isfcalc,& end do read_loop end do read_subset call closbf(lnbufr) + close(lnbufr) deallocate(data1b8) num_obs = iob-1 diff --git a/src/gsi/read_satmar.f90 b/src/gsi/read_satmar.f90 index 3a7f349c2b..93a7780b1b 100644 --- a/src/gsi/read_satmar.f90 +++ b/src/gsi/read_satmar.f90 @@ -217,7 +217,7 @@ subroutine read_satmar (nread, ndata, nodata, & endif ! ! *#* Main - Start *#*! - call closbf(lun11) +! call closbf(lun11) open(lun11,file=trim(infile),action='read',form='unformatted', iostat=ierr) if (ierr/=0) then print*, myname,' : ERROR : File ', trim(infile),' not existing. ' @@ -234,6 +234,7 @@ subroutine read_satmar (nread, ndata, nodata, & end do end do call closbf(lun11) + close(lun11) ! ! Allocate Arrays for all the data allocate (data_all (nreal, cnt),isort(cnt)) @@ -491,7 +492,7 @@ subroutine read_satmar (nread, ndata, nodata, & deallocate(data_out) if (ndata == 0) then - write(6,*)myname,': closbf(',lun11,')' + write(6,*)myname,': closbf(',lun11,') no data' endif close(lun11) ! diff --git a/src/gsi/read_satwnd.f90 b/src/gsi/read_satwnd.f90 index 1bc4f1b8dc..2d22aeb3ab 100644 --- a/src/gsi/read_satwnd.f90 +++ b/src/gsi/read_satwnd.f90 @@ -304,7 +304,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis nmsg=0 nrep=0 ntb =0 - call closbf(lunin) +! call closbf(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -589,6 +589,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis endif call closbf(lunin) + close(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -1544,8 +1545,6 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis enddo loop_readsb ! End of bufr read loop 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) @@ -1560,6 +1559,8 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis enddo loop_convinfo! loops over convinfo entry matches deallocate(lmsg,tab,nrep) +! Close unit to bufr file + call closbf(lunin) ! Write header record and data to output file for further processing @@ -1600,10 +1601,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_SATWND: ',& 'nvtest,vdisterrmax=',ntest,vdisterrmax - - if (ndata == 0) then - call closbf(lunin) write(6,*)'READ_SATWND: closbf(',lunin,')' endif diff --git a/src/gsi/read_seviri.f90 b/src/gsi/read_seviri.f90 index 9b1eb63657..853b69620f 100644 --- a/src/gsi/read_seviri.f90 +++ b/src/gsi/read_seviri.f90 @@ -190,7 +190,7 @@ subroutine read_seviri(mype,val_sev,ithin,rmesh,jsatid,& endif ! Open bufr file. - call closbf(lnbufr) +! call closbf(lnbufr) open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) @@ -247,6 +247,7 @@ subroutine read_seviri(mype,val_sev,ithin,rmesh,jsatid,& ! Reopen unit to bufr file call closbf(lnbufr) + close(lnbufr) if(jsatid == 'm08') kidsat = 55 if(jsatid == 'm09') kidsat = 56 if(jsatid == 'm10') kidsat = 57 @@ -526,6 +527,7 @@ subroutine read_seviri(mype,val_sev,ithin,rmesh,jsatid,& if(allocated(rd_tdiffs)) deallocate(rd_tdiffs) enddo read_msg call closbf(lnbufr) + close(lnbufr) call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& nele,itxmax,nread,ndata,data_all,score_crit,nrec) @@ -621,6 +623,7 @@ subroutine read_subset_nnsb subset_num(nnmsg)=nnsb enddo read_msg2 call closbf(lnbufr) + close(lnbufr) if (npe_sub > 1 ) then call mpi_allreduce(subset_num, subset_nnsb, nnmsg,mpi_itype,mpi_sum,mpi_comm_sub,ierror) else diff --git a/src/gsi/read_sfcwnd.f90 b/src/gsi/read_sfcwnd.f90 index cabe69b84f..ac1174fef3 100644 --- a/src/gsi/read_sfcwnd.f90 +++ b/src/gsi/read_sfcwnd.f90 @@ -231,7 +231,7 @@ subroutine read_sfcwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis !! go through the satedump to find out how many subset to process !** Open and read data from bufr data file - call closbf(lunin) +! call closbf(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -379,6 +379,7 @@ subroutine read_sfcwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis endif call closbf(lunin) + close(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -724,7 +725,6 @@ subroutine read_sfcwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis 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) @@ -733,6 +733,7 @@ subroutine read_sfcwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! Normal exit enddo loop_convinfo! loops over convinfo entry matches + call closbf(lunin) deallocate(lmsg,nrep,tab) @@ -772,8 +773,7 @@ subroutine read_sfcwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis 'nvtest,vdisterrmax=',ntest,vdisterrmax if (ndata == 0) then - call closbf(lunin) - write(6,*)'READ_SFCWND: closbf(',lunin,')' + write(6,*)'READ_SFCWND: closbf(',lunin,') no data' endif write(6,*) 'READ_SFCWND,nread,ndata,nreal,nodata=',nread,ndata,nreal,nodata diff --git a/src/gsi/read_ssmi.f90 b/src/gsi/read_ssmi.f90 index 82c7c39ae5..3e47ee79b5 100755 --- a/src/gsi/read_ssmi.f90 +++ b/src/gsi/read_ssmi.f90 @@ -511,6 +511,7 @@ subroutine read_ssmi(mype,val_ssmi,ithin,rmesh,jsatid,gstime,& end do read_loop end do read_subset call closbf(lnbufr) + close(lnbufr) ! If multiple tasks read input bufr file, allow each tasks to write out ! information it retained and then let single task merge files together diff --git a/src/gsi/read_ssmis.f90 b/src/gsi/read_ssmis.f90 index faaec0cece..3e0da712f6 100755 --- a/src/gsi/read_ssmis.f90 +++ b/src/gsi/read_ssmis.f90 @@ -390,7 +390,7 @@ subroutine read_ssmis(mype,val_ssmis,ithin,isfcalc,rmesh,jsatid,gstime,& ! Read in data from bufr into arrays first ! Open unit to satellite bufr file iobs=1 - call closbf(lnbufr) +! call closbf(lnbufr) open(lnbufr,file=trim(infile),form='unformatted',status='old',err=500) call openbf(lnbufr,'IN',lnbufr) call datelen(10) @@ -516,6 +516,7 @@ subroutine read_ssmis(mype,val_ssmis,ithin,isfcalc,rmesh,jsatid,gstime,& end do read_loop end do read_subset call closbf(lnbufr) + close(lnbufr) num_obs = iobs-1 diff --git a/src/gsi/read_wcpbufr.f90 b/src/gsi/read_wcpbufr.f90 index 7d94e7bd1a..560d58134e 100644 --- a/src/gsi/read_wcpbufr.f90 +++ b/src/gsi/read_wcpbufr.f90 @@ -208,7 +208,7 @@ subroutine read_wcpbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& irec = 0 ! Open, then read date from bufr data - call closbf(lunin) +! call closbf(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -342,6 +342,7 @@ subroutine read_wcpbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& call closbf(lunin) + close(lunin) open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -664,8 +665,6 @@ subroutine read_wcpbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& end do loop_readsb ! End of bufr read loop enddo loop_msg -! Close unit to bufr file - call closbf(lunin) ! Deallocate arrays used for thinning data if (.not.use_all) then @@ -676,6 +675,9 @@ subroutine read_wcpbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! Normal exit enddo loop_convinfo! loops over convinfo entry matches +! Close unit to bufr file + call closbf(lunin) + close(lunin) deallocate(lmsg,tab,nrep) ! Write header record and data to output file for further processing @@ -713,8 +715,7 @@ subroutine read_wcpbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& 'nvtest,vdisterrmax=',ntest,vdisterrmax if (ndata == 0) then - call closbf(lunin) - write(6,*)'READ_WCPBUFR: closbf(',lunin,')' + write(6,*)'READ_WCPBUFR: closbf(',lunin,') no data' endif close(lunin) diff --git a/src/gsi/setupps.f90 b/src/gsi/setupps.f90 index 3851b8d10e..65e6fe999c 100644 --- a/src/gsi/setupps.f90 +++ b/src/gsi/setupps.f90 @@ -133,6 +133,7 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa use m_dtime, only: dtime_setup, dtime_check + use hdraobmod, only: nhdps,hdpslist use gsi_bundlemod, only : gsi_bundlegetpointer use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle use sparsearr, only: sparr2, new, size, writearray, fullarray @@ -180,9 +181,9 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa integer(i_kind) ier,ilon,ilat,ipres,ihgt,itemp,id,itime,ikx,iqc,iptrb,ijb integer(i_kind) ier2,iuse,ilate,ilone,istnelv,idomsfc,izz,iprvd,isprvd integer(i_kind) ikxx,nn,ibin,ioff,ioff0 - integer(i_kind) i,nchar,nreal,ii,jj,k,l,mm1 + integer(i_kind) i,j,nchar,nreal,ii,jj,k,l,mm1 integer(i_kind) itype,isubtype - integer(i_kind) ibb,ikk + integer(i_kind) ibb,ikk,idddd logical,dimension(nobs):: luse,muse integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID @@ -268,6 +269,25 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa do i=1,nobs muse(i)=nint(data(iuse,i)) <= jiter end do +! If HD raobs available move prepbufr version to monitor + if(nhdps > 0)then + do i=1,nobs + ikx=nint(data(ikxx,i)) + itype=ictype(ikx) + if(itype == 120) then + rstation_id = data(id,i) + read(station_id,'(i5,3x)',err=1200) idddd + stn_loop:do j=1,nhdps + if(idddd == hdpslist(j))then + data(iuse,i)=108. + muse(i) = .false. + exit stn_loop + end if + end do stn_loop + end if +1200 continue + end do + end if hr_offset=min_offset/60.0_r_kind ! Check for duplicate observations at same location diff --git a/src/gsi/setupq.f90 b/src/gsi/setupq.f90 index c445bba29f..d9ce9ec14c 100644 --- a/src/gsi/setupq.f90 +++ b/src/gsi/setupq.f90 @@ -169,6 +169,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle use sparsearr, only: sparr2, new, size, writearray, fullarray use state_vectors, only: svars3d, levels + use hdraobmod, only: nhdq,hdqlist ! The following variables are the coefficients that describe the ! linear regression fits that are used to define the dynamic @@ -231,15 +232,15 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav real(r_single),allocatable,dimension(:,:)::rdiagbufp - integer(i_kind) i,nchar,nreal,ii,l,jj,mm1,itemp,iip + integer(i_kind) i,j,nchar,nreal,ii,l,jj,mm1,itemp,iip integer(i_kind) jsig,itype,k,nn,ikxx,iptrb,ibin,ioff,ioff0,icat,ijb integer(i_kind) ier,ilon,ilat,ipres,iqob,id,itime,ikx,iqmax,iqc integer(i_kind) ier2,iuse,ilate,ilone,istnelv,iobshgt,izz,iprvd,isprvd integer(i_kind) idomsfc,iderivative - integer(i_kind) ibb,ikk + integer(i_kind) ibb,ikk,idddd real(r_kind) :: delz type(sparr2) :: dhx_dx - integer(i_kind) :: iz, q_ind, nind, nnz + integer(i_kind) :: iz, q_ind, nind, nnz,iprev_station character(8) station_id character(8),allocatable,dimension(:):: cdiagbuf,cdiagbufp @@ -317,6 +318,34 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav do i=1,nobs muse(i)=nint(data(iuse,i)) <= jiter end do +! If HD raobs available move prepbufr version to monitor + nhdq=0 + if(nhdq > 0)then + iprev_station=0 + do i=1,nobs + ikx=nint(data(ikxx,i)) + itype=ictype(ikx) + if(itype == 120) then + rstation_id = data(id,i) + read(station_id,'(i5,3x)',err=1200) idddd + if(idddd == iprev_station)then + data(iuse,i)=108. + muse(i) = .false. + else + stn_loop:do j=1,nhdq + if(idddd == hdqlist(j))then + iprev_station=idddd + data(iuse,i)=108. + muse(i) = .false. +! write(6,*) ' in setupq ',idddd + exit stn_loop + end if + end do stn_loop + end if + end if +1200 continue + end do + end if var_jb=zero @@ -397,12 +426,12 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Prepare specific humidity data call dtime_setup() do i=1,nobs + ikx=nint(data(ikxx,i)) + itype=ictype(ikx) dtime=data(itime,i) call dtime_check(dtime, in_curbin, in_anybin) if(.not.in_anybin) cycle - ikx=nint(data(ikxx,i)) - itype=ictype(ikx) ! Flag static conditions to create PBL_pseudo_surfobsq obs. l_pbl_pseudo_itype = l_pbl_pseudo_surfobsq .and. & @@ -415,7 +444,6 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav dpres=data(ipres,i) rmaxerr=data(iqmax,i) - rstation_id = data(id,i) error=data(ier2,i) prest=r10*exp(dpres) ! in mb var_jb=data(ijb,i) @@ -1059,6 +1087,7 @@ subroutine init_netcdf_diag_ end subroutine init_netcdf_diag_ subroutine contents_binary_diag_(odiag) type(obs_diag),pointer,intent(in):: odiag + cdiagbuf(ii) = station_id ! station id rdiagbuf(1,ii) = ictype(ikx) ! observation type diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index 06d50497c3..71b6867506 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -79,6 +79,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use gsi_bundlemod, only : gsi_bundlegetpointer use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle use buddycheck_mod, only: buddy_check_t + use hdraobmod, only: nhdt,hdtlist use sparsearr, only: sparr2, new, size, writearray, fullarray @@ -285,13 +286,13 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav integer(i_kind) ier2,iuse,ilate,ilone,ikxx,istnelv,iobshgt,izz,iprvd,isprvd integer(i_kind) regime integer(i_kind) idomsfc,iskint,iff10,isfcr - integer(i_kind) ibb,ikk + integer(i_kind) ibb,ikk,idddd integer(i_kind),dimension(nobs):: buddyuse type(sparr2) :: dhx_dx - integer(i_kind) :: iz, t_ind, nind, nnz + integer(i_kind) :: iz, t_ind, nind, nnz, iprev_station character(8) station_id character(8),allocatable,dimension(:):: cdiagbuf,cdiagbufp character(8),allocatable,dimension(:):: cprvstg,csprvstg @@ -391,6 +392,33 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav do i=1,nobs muse(i)=nint(data(iuse,i)) <= jiter end do +! If HD raobs available move prepbufr version to monitor + if(nhdt > 0)then + iprev_station=0 + do i=1,nobs + ikx=nint(data(ikxx,i)) + itype=ictype(ikx) + if(itype == 120) then + rstation_id = data(id,i) + read(station_id,'(i5,3x)',err=1200) idddd + if(idddd == iprev_station)then + data(iuse,i)=108. + muse(i) = .false. + else + stn_loop:do j=1,nhdt + if(idddd == hdtlist(j))then + iprev_station=idddd + data(iuse,i)=108. + muse(i) = .false. +! write(6,*) ' in setupt ',idddd + exit stn_loop + end if + end do stn_loop + end if + end if +1200 continue + end do + end if var_jb=zero ! handle multiple reported data at a station diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90 index d14f2c0fd8..f73cd416ab 100644 --- a/src/gsi/setupw.f90 +++ b/src/gsi/setupw.f90 @@ -69,10 +69,12 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use m_dtime, only: dtime_setup, dtime_check use gsi_bundlemod, only : gsi_bundlegetpointer + use hdraobmod, only: nhduv,hduvlist use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle use sparsearr, only: sparr2, new, size, writearray, fullarray use aux2dvarflds, only: rtma_comp_fact10 + ! The following variables are the coefficients that describe the ! linear regression fits that are used to define the dynamic ! observation error (DOE) specifications for all reconnissance @@ -285,9 +287,9 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav integer(i_kind) ihgt,ier2,iuse,ilate,ilone integer(i_kind) izz,iprvd,isprvd integer(i_kind) idomsfc,isfcr,iskint,iff10 - integer(i_kind) ibb,ikk,ihil + integer(i_kind) ibb,ikk,ihil,idddd - integer(i_kind) num_bad_ikx + integer(i_kind) num_bad_ikx,iprev_station character(8) station_id character(8),allocatable,dimension(:):: cdiagbuf @@ -409,9 +411,42 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if (netcdf_diag) call init_netcdf_diag_ end if + num_bad_ikx=0 do i=1,nobs muse(i)=nint(data(iuse,i)) <= jiter + ikx=nint(data(ikxx,i)) + if(ikx < 1 .or. ikx > nconvtype) then + num_bad_ikx=num_bad_ikx+1 + if(num_bad_ikx<=10) write(6,*)' in setupw ',ikx,i,nconvtype,mype + end if end do +! If HD raobs available move prepbufr version to monitor + if(nhduv > 0)then + iprev_station=0 + do i=1,nobs + ikx=nint(data(ikxx,i)) + itype=ictype(ikx) + if(itype == 220 .or. itype == 221) then + rstation_id = data(id,i) + read(station_id,'(i5,3x)',err=1200) idddd + if(idddd == iprev_station)then + data(iuse,i)=108. + muse(i) = .false. + else + stn_loop:do j=1,nhduv + if(idddd == hduvlist(j))then + iprev_station=idddd + data(iuse,i)=108. + muse(i) = .false. +! write(6,*) ' in setupw ',idddd,itype + exit stn_loop + end if + end do stn_loop + end if + end if +1200 continue + end do + end if ! handle multiple-report observations at a station hr_offset=min_offset/60.0_r_kind @@ -468,7 +503,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav var_jb=data(ijb,i) if(ikx < 1 .or. ikx > nconvtype) then num_bad_ikx=num_bad_ikx+1 - if(num_bad_ikx<=10) write(6,*)' in setupw, bad ikx, ikx,i,nconvtype=',ikx,i,nconvtype + if(num_bad_ikx<=10) write(6,*)' in setupw, bad ikx, ikx,i,nconvtype=',ikx,i,nconvtype,mype cycle end if isli = data(idomsfc,i) @@ -552,6 +587,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav z_height = .false. if ((itype>=221 .and. itype <= 229) .and. (data(ihgt,i)=223 .and. itype<=228) .or. sfc_data) .and. .not.twodvar_regional) then + if (((itype>=223 .and. itype<=228) .or. itype == 218 .or. sfc_data) .and. .not.twodvar_regional) then ! Convert geopotential height at layer midpoints to geometric ! height using equations (17, 20, 23) in MJ Mahoney's note ! "A discussion of various measures of altitude" (2001). diff --git a/src/gsi/sst_retrieval.f90 b/src/gsi/sst_retrieval.f90 index bf530b29cb..97f1073a50 100644 --- a/src/gsi/sst_retrieval.f90 +++ b/src/gsi/sst_retrieval.f90 @@ -468,6 +468,7 @@ subroutine finish_sst_retrieval implicit none call closbf(lnbufr) + close(lnbufr) return end subroutine finish_sst_retrieval diff --git a/ush/build_all_cmake.sh b/ush/build_all_cmake.sh index 0205eeccb4..02f46cf7b5 100755 --- a/ush/build_all_cmake.sh +++ b/ush/build_all_cmake.sh @@ -89,6 +89,6 @@ else cmake .. fi -make -j 8 +make -j 8 exit