From 6e0da830f0a2a4af2fe88735810891fc32e6a344 Mon Sep 17 00:00:00 2001 From: JingCheng-NOAA Date: Mon, 28 Aug 2023 18:57:41 +0000 Subject: [PATCH] Merge HAFS GSI changes with latest GSI dev branch --- src/gsi/gsi_obOperTypeManager.F90 | 8 +- src/gsi/m_extOzone.F90 | 319 ++++++++++++++++++++++++++++-- src/gsi/read_bufrtovs.f90 | 2 +- src/gsi/read_obs.F90 | 10 +- src/gsi/read_satwnd.f90 | 92 +++++++-- src/gsi/read_ssmi.f90 | 6 +- src/gsi/setupoz.f90 | 156 +++++++++++---- 7 files changed, 507 insertions(+), 86 deletions(-) diff --git a/src/gsi/gsi_obOperTypeManager.F90 b/src/gsi/gsi_obOperTypeManager.F90 index ea306953c4..5df899825a 100644 --- a/src/gsi/gsi_obOperTypeManager.F90 +++ b/src/gsi/gsi_obOperTypeManager.F90 @@ -276,6 +276,9 @@ function dtype2index_(dtype) result(index_) case("ompstc8"); index_= iobOper_oz case("ompsnp" ); index_= iobOper_oz case("ompsnm" ); index_= iobOper_oz + case("omieff" ); index_= iobOper_oz + case("tomseff" ); index_= iobOper_oz + case("ompsnmeff"); index_= iobOper_oz case("o3l" ,"[o3loper]" ); index_= iobOper_o3l case("o3lev" ); index_= iobOper_o3l @@ -283,11 +286,10 @@ function dtype2index_(dtype) result(index_) case("mls22" ); index_= iobOper_o3l case("mls30" ); index_= iobOper_o3l case("mls55" ); index_= iobOper_o3l - case("omieff" ); index_= iobOper_o3l - case("tomseff" ); index_= iobOper_o3l + case("ompslp" ); index_= iobOper_o3l case("ompslpuv" ); index_= iobOper_o3l case("ompslpvis"); index_= iobOper_o3l - case("ompslp" ); index_= iobOper_o3l + case("ompslpnc" ); index_= iobOper_o3l case("gpsbend","[gpsbendoper]"); index_= iobOper_gpsbend case("gps_bnd"); index_= iobOper_gpsbend diff --git a/src/gsi/m_extOzone.F90 b/src/gsi/m_extOzone.F90 index 3d4b6783c1..a28209292f 100644 --- a/src/gsi/m_extOzone.F90 +++ b/src/gsi/m_extOzone.F90 @@ -158,32 +158,36 @@ function is_extOzone_(dfile,dtype,dplat,class) is_extOzone_= & ifile_==iBUFR .and. dtype == 'o3lev' .or. & ifile_==iNC .and. dtype == 'mls55' .or. & + ifile_==iNC .and. dtype == 'ompslpnc' .or. & ifile_==iNC .and. dtype == 'ompslpvis' .or. & ifile_==iNC .and. dtype == 'ompslpuv' .or. & - ifile_==iNC .and. dtype == 'ompslp' .or. & ifile_==iNC .and. dtype == 'lims' .or. & ifile_==iNC .and. dtype == 'uarsmls' .or. & ifile_==iNC .and. dtype == 'mipas' .or. & ifile_==iNC .and. dtype == 'omieff' .or. & + ifile_==iNC .and. dtype == 'ompsnmeff' .or. & + ifile_==iNC .and. dtype == 'ompsnpnc' .or. & ifile_==iNC .and. dtype == 'tomseff' case(iLEVEL) is_extOzone_= & ifile_==iBUFR .and. dtype == 'o3lev' .or. & ifile_==iNC .and. dtype == 'mls55' .or. & + ifile_==iNC .and. dtype == 'ompslpnc' .or. & ifile_==iNC .and. dtype == 'ompslpvis' .or. & ifile_==iNC .and. dtype == 'ompslpuv' .or. & - ifile_==iNC .and. dtype == 'ompslp' .or. & ifile_==iNC .and. dtype == 'lims' .or. & ifile_==iNC .and. dtype == 'uarsmls' .or. & ifile_==iNC .and. dtype == 'mipas' case(iLAYER) - is_extOzone_= .false. + is_extOzone_= & + ifile_==iNC .and. dtype == 'ompsnpnc' case(iTOTAL) is_extOzone_= & - ifile_==iNC .and. dtype == 'omieff' .or. & + ifile_==iNC .and. dtype == 'omieff' .or. & + ifile_==iNC .and. dtype == 'ompsnmeff' .or. & ifile_==iNC .and. dtype == 'tomseff' case default @@ -332,7 +336,7 @@ subroutine read_(dfile,dtype,dplat,dsis, & ! intent(in), keys for type mana endif select case(dtype) - case('omieff','tomseff') ! layer-ozone or total-ozone types + case('omieff','tomseff','ompsnmeff') ! layer-ozone or total-ozone types select case(dfile_format(dfile)) case('nc') call oztot_ncInquire_(nreal,nchan,ilat,ilon, & @@ -381,7 +385,7 @@ subroutine read_(dfile,dtype,dplat,dsis, & ! intent(in), keys for type mana jsatid, gstime,twind) end select - case('mls55','ompslpvis','ompslpuv','ompslp','lims','uarsmls','mipas') + case('mls55','ompslpnc','ompslpvis','ompslpuv','lims','uarsmls','mipas') select case(dfile_format(dfile)) case('nc') call ozlev_ncInquire_( nreal,nchan,ilat,ilon,maxobs) @@ -393,6 +397,17 @@ subroutine read_(dfile,dtype,dplat,dsis, & ! intent(in), keys for type mana end select + case('ompsnpnc') + select case(dfile_format(dfile)) + case('nc') + call ozlay_ncInquire_( nreal,nchan,ilat,ilon,maxobs) + allocate(p_out(nreal+nchan,maxobs)) + p_out(:,:)=RMISS + + call ozlay_ncRead_(dfile,dtype, p_out,nread,npuse,nouse, gstime,twind) + + end select + end select if(nouse<0 .or. .not.associated(p_out)) then @@ -706,7 +721,7 @@ subroutine oztot_ncread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, & ! Apply data screening based on quality flags ! Bit 10 (from the left) in TOQF represents row anomaly. All 17 bits in toqf is ! supposed to converted into array elements of binary(:), either for "tomseff" or -! "omieff". +! "omieff" or "ompsnmeff". binary(:)=0 select case(dtype) case('omieff') @@ -731,6 +746,9 @@ subroutine oztot_ncread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, & ! 0 - good data, 1 - good data with SZA > 84 deg if (toqf /= 0) cycle recloop + case('ompsnmeff') + !! data in NetCDF are prescreened + case default end select @@ -764,10 +782,10 @@ subroutine oztot_ncread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, & ozout(7,itx)=real(toqf) ! - total ozone quality code (not used) ozout(8,itx)=real(sza) ! solar zenith angle ozout(9,itx)=binary(10) ! row anomaly flag, is actually fixed to 0 - ozout(10,itx)=0. ! - cloud amount (not used) - ozout(11,itx)=0. ! - vzan (not used) - ozout(12,itx)=0. ! - aerosol index (not used) - ozout(13,itx)=0. ! - ascending/descending (not used) + ozout(10,itx)=zero ! - cloud amount (not used) + ozout(11,itx)=zero ! - vzan (not used) + ozout(12,itx)=zero ! - aerosol index (not used) + ozout(13,itx)=zero ! - ascending/descending (not used) ozout(14,itx)=real(fovn) ! scan position ! "(not used)" flags above imply that they ! are not used in setupozlay(). @@ -1421,13 +1439,286 @@ 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 -end subroutine ozlev_bufrread_ + end subroutine ozlev_bufrread_ + + subroutine ozlay_ncInquire_( nreal,nchan,ilat,ilon, maxrec) + implicit none + + integer(kind=i_kind), intent(out):: nreal ! number of real parameters per record + integer(kind=i_kind), intent(out):: nchan ! number of channels or levels per record + integer(kind=i_kind), intent(out):: ilat ! index to latitude in nreal parameters. + integer(kind=i_kind), intent(out):: ilon ! index to longitude in nreal parameters. + + integer(kind=i_kind), intent(out):: maxrec ! extimated input record count + + character(len=*), parameter:: myname_=myname//'::ozlay_ncInquire_' + + ! Configure the record, they are not (dfile,dtype,dplat) dependent in this case. + nreal = 9 + nchan = 22 + ilat=4 + ilon=3 + + maxrec = MAXOBS_ + end subroutine ozlay_ncInquire_ + + !.................................................................................. + subroutine ozlay_ncread_(dfile,dtype,ozout,nmrecs,ndata,nodata, gstime,twind) + !.................................................................................. + use netcdf, only: nf90_open + use netcdf, only: nf90_nowrite + use netcdf, only: nf90_noerr + use netcdf, only: nf90_inq_dimid + use netcdf, only: nf90_inquire_dimension + use netcdf, only: nf90_inq_varid + use netcdf, only: nf90_get_var + use netcdf, only: nf90_close + + use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons + use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar + + use constants, only: deg2rad,zero,rad2deg,one_tenth,r60inv + use ozinfo, only: jpch_oz,nusis_oz,iuse_oz + use mpeu_util, only: perr,die + ! use mpeu_util, only: mprefix,stdout + + implicit none + character(len=*), intent(in):: dfile ! obs_input filename + character(len=*), intent(in):: dtype ! obs_input dtype + real (kind=r_kind), dimension(:,:), intent(out):: ozout + integer(kind=i_kind), intent(out):: nmrecs ! count of records read + integer(kind=i_kind), intent(out):: ndata ! count of processed + integer(kind=i_kind), intent(out):: nodata ! count of retained + real (kind=r_kind), intent(in):: gstime ! analysis time (minute) from reference date + real (kind=r_kind), intent(in):: twind ! input group time window (hour) + + character(len=*), parameter:: myname_=myname//'::ozlay_ncRead_' + + integer(kind=i_kind):: ier,nprofs,levs,ikx,i,k0,ilev,iprof + integer(kind=i_kind):: nrecDimId,lonVarID,latVarID,yyVarID,mmVarID,levsDimId + integer(kind=i_kind):: szaVarID,ozoneVarID,nmind,k + integer(kind=i_kind):: ddVarID,hhVarID,minVarID,ssVarID,maxobs,ncid + real (kind=r_kind):: dlon,dlon_earth,dlon_earth_deg + real (kind=r_kind):: dlat,dlat_earth,dlat_earth_deg + real (kind=r_kind):: slons0,slats0 + real (kind=r_kind):: tdiff,sstime,t4dv,rsat + integer(kind=i_kind):: idate5(5) + integer(kind=i_kind),allocatable,dimension(:):: ipos + real(r_kind),allocatable,dimension(:):: poz + + integer(kind=i_kind), allocatable :: iya(:),ima(:),idda(:),ihha(:),imina(:),iseca(:) + real (kind=r_kind), allocatable :: slatsa(:),slonsa(:),ozone(:,:),sza(:) + real(r_kind) totoz + + logical:: outside + logical:: first + real(r_kind),parameter:: badoz = 10000.0_r_kind + + maxobs=size(ozout,2) + rsat=999._r_kind + ndata = 0 + nmrecs=0 + nodata=-1 + + ! Open file and read dimensions + call check(nf90_open(trim(dfile),nf90_nowrite,ncid),stat=ier) + + ! ignore if the file is not actually present. + if(ier/=nf90_noerr) then + nodata = 0 + return + endif + + ! Get dimensions from the input file + call check(nf90_inq_dimid(ncid, "nrec", nrecDimId),stat=ier) + + ! ignore if error + if(ier/=nf90_noerr) then + nodata = 0 + call check(nf90_close(ncid),stat=ier) + return + endif + + ! Get dimensions from the input file: # of profiles and # of levels + nprofs=0 + call check(nf90_inquire_dimension(ncid, nrecDimId, len = nprofs),stat=ier) + ! ignore if error + if(ier/=nf90_noerr) then + call check(nf90_close(ncid),stat=ier) + return + endif + + if(nprofs==0) then + nodata=0 + call check(nf90_close(ncid),stat=ier) + return + endif + + ! Continue the input + call check(nf90_inq_dimid(ncid, "nlevs", levsDimId)) + call check(nf90_inquire_dimension(ncid, levsDimId, len = levs)) + !!!!! if (levs /= nchan) + + allocate(ipos(levs)) + ipos=999 + ikx = 0 + first=.false. + do i=1,jpch_oz + if( (.not. first) .and. index(nusis_oz(i), trim(dtype))/=0) then + k0=i + first=.true. + end if + if(first.and.index(nusis_oz(i),trim(dtype))/=0) then + ikx=ikx+1 + ipos(ikx)=k0+ikx-1 + end if + end do + + if (ikx/=levs+1) call die(myname_//': inconsistent levs for '//dtype) + + ! Allocate space and read data + allocate(iya(nprofs),ima(nprofs),idda(nprofs),ihha(nprofs),imina(nprofs), & + iseca(nprofs),slatsa(nprofs),slonsa(nprofs),sza(nprofs),ozone(levs,nprofs)) + allocate (poz(levs+1)) + + call check(nf90_inq_varid(ncid, "lon", lonVarId)) + call check(nf90_get_var(ncid, lonVarId, slonsa)) + + call check(nf90_inq_varid(ncid, "lat", latVarId)) + call check(nf90_get_var(ncid, latVarId, slatsa)) + + call check(nf90_inq_varid(ncid, "yy", yyVarId)) + call check(nf90_get_var(ncid, yyVarId, iya)) + + call check(nf90_inq_varid(ncid, "mm", mmVarId)) + call check(nf90_get_var(ncid, mmVarId, ima)) + + call check(nf90_inq_varid(ncid, "dd", ddVarId)) + call check(nf90_get_var(ncid, ddVarId, idda)) + + call check(nf90_inq_varid(ncid, "hh", hhVarId)) + call check(nf90_get_var(ncid, hhVarId, ihha)) + + call check(nf90_inq_varid(ncid, "min", minVarId)) + call check(nf90_get_var(ncid, minVarId, imina)) + + call check(nf90_inq_varid(ncid, "ss", ssVarId)) + call check(nf90_get_var(ncid, ssVarId, iseca)) + + call check(nf90_inq_varid(ncid, "sza", szaVarId)) + call check(nf90_get_var(ncid, szaVarId, sza)) + + call check(nf90_inq_varid(ncid, "ozone", ozoneVarId)) + call check(nf90_get_var(ncid, ozoneVarId, ozone)) + + ! close the data file + call check(nf90_close(ncid)) + + ! 'Unpack' the data + nmrecs = 0 + nodata = 0 + read_loop1: do iprof = 1, nprofs + do ilev = 1, levs + if (ozone(ilev, iprof) .lt. -900.0_r_kind) cycle ! undefined + end do +!!$ if (iuse_oz(ipos(ilev)) < 0) then +!!$ usage = 10000._r_kind +!!$ else +!!$ usage = zero +!!$ endif + nmrecs=nmrecs+levs+1 + + ! convert observation location to radians + slons0=slonsa(iprof) + slats0=slatsa(iprof) + if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle + if(slons0< zero) slons0=slons0+r360 + if(slons0==r360) slons0=zero + dlat_earth_deg = slats0 + dlon_earth_deg = slons0 + dlat_earth = slats0 * deg2rad + dlon_earth = slons0 * deg2rad + + if(regional)then + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) + if(outside) cycle + else + dlat = dlat_earth + dlon = dlon_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + + idate5(1) = iya(iprof) !year + idate5(2) = ima(iprof) !month + idate5(3) = idda(iprof) !day + idate5(4) = ihha(iprof) !hour + idate5(5) = imina(iprof) !minute + call w3fs21(idate5,nmind) + t4dv=real((nmind-iwinbgn),r_kind)*r60inv + if (l4dvar.or.l4densvar) then + if (t4dvwinlen) then + write(6,*)'read_ozone: ', dtype,' obs time idate5=',idate5,', t4dv=',& + t4dv,' is outside time window, sstime=',sstime*r60inv + cycle + end if + else + sstime=real(nmind,r_kind) + tdiff=(sstime-gstime)*r60inv + if(abs(tdiff) > twind)then + write(6,*)'read_ozone: ',dtype,' obs time idate5=',idate5,', tdiff=',& + tdiff,' is outside time window=',twind + cycle + end if + end if + + !! Compute total ozone + totoz=zero + do k=1,levs + poz(k) = ozone(k,iprof) + totoz=totoz+ozone(k,iprof) + end do + poz(levs+1) = totoz + + !Check ozone layer values. If any layer value is bad, toss entire profile + do k=1,levs + if (poz(k)>badoz) cycle read_loop1 + end do + + ! Write ozone record to output file + ndata=min(ndata+1,maxobs) + if(ndata<=maxobs) then + nodata=nodata+levs+1 + ozout(1,ndata)=rsat + ozout(2,ndata)=t4dv + ozout(3,ndata)=dlon ! grid relative longitude + ozout(4,ndata)=dlat ! grid relative latitude + ozout(5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) + ozout(6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) + ozout(7,ndata)=zero ! total ozone error flag + ozout(8,ndata)=zero ! profile ozone error flag + ozout(9,ndata)=sza(iprof) ! solar zenith angle + do k=1,levs+1 + ozout(k+9,ndata)=poz(k) + end do + end if + end do read_loop1 + + deallocate(iya,ima,idda,ihha,imina,iseca,slatsa,slonsa, ozone, poz,sza) + deallocate(ipos) + if (ndata > maxobs) then + call perr('read_ozone','Number of layer obs reached maxobs = ', maxobs) + call perr(myname_,'Number of layer obs reached maxobs = ', maxobs) + call perr(myname_,' ndata = ', ndata) + call perr(myname_,' nodata = ', nodata) + call die(myname_) + endif + + end subroutine ozlay_ncread_ !.................................................................................. subroutine check(status,stat) diff --git a/src/gsi/read_bufrtovs.f90 b/src/gsi/read_bufrtovs.f90 index 0c954c7c1d..a819acd2c3 100644 --- a/src/gsi/read_bufrtovs.f90 +++ b/src/gsi/read_bufrtovs.f90 @@ -683,7 +683,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& if (llll > 1) then sacv = nint(bfr1bhdr(14)) if (sacv > spc_coeff_versions) then - write(6,*) 'READ_BUFRTOVS WARNING sacv greater than spc_coeff_versions' + write(6,*) 'READ_BUFRTOVS WARNING sacv greater than spc_coeff_versions',' ',jsatid,' ',obstype end if else ! normal feed doesn't have antenna correction, so set sacv to 0 sacv = 0 diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index c9eeb8a869..e564e323fc 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -892,6 +892,7 @@ subroutine read_obs(ndata,mype) if(obstype == 'mls20' ) nmls_type=nmls_type+1 if(obstype == 'mls22' ) nmls_type=nmls_type+1 if(obstype == 'mls30' ) nmls_type=nmls_type+1 + if(obstype == 'mls55' ) nmls_type=nmls_type+1 if(nmls_type>1) then write(6,*) '******ERROR***********: there is more than one MLS data type, not allowed, please check' call stop2(339) @@ -935,6 +936,7 @@ subroutine read_obs(ndata,mype) .or. obstype == 'ompsnp' & .or. obstype == 'gome' & .or. index(obstype, 'omps') /= 0 & + .or. index(obstype, 'omi' ) /= 0 & .or. mls & ) then ditype(i) = 'ozone' @@ -1081,8 +1083,12 @@ subroutine read_obs(ndata,mype) if (ii>npem1) ii=0 if(mype==ii)then call gsi_inquire(lenbytes,lexist,trim(dfile(i)),mype) - call read_obs_check (lexist,trim(dfile(i)),dplat(i),dtype(i),minuse,read_rec1(i)) - + !call read_obs_check (lexist,trim(dfile(i)),dplat(i),dtype(i),minuse,read_rec1(i)) + if (is_extOzone(dfile(i),obstype,dplat(i))) then + print*,'reading ',trim(dfile(i)),' ',obstype,' ',trim(dplat(i)),lexist,lenbytes + else + call read_obs_check(lexist,trim(dfile(i)),dplat(i),dtype(i),minuse,read_rec1(i)) + endif ! If no data set starting record to be 999999. Note if this is not large ! enough code should still work - just does a bit more work. diff --git a/src/gsi/read_satwnd.f90 b/src/gsi/read_satwnd.f90 index d50da90a83..dc1bfbb34e 100644 --- a/src/gsi/read_satwnd.f90 +++ b/src/gsi/read_satwnd.f90 @@ -158,11 +158,11 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis real(r_kind),parameter:: r1200= 1200.0_r_kind real(r_kind),parameter:: r10000= 10000.0_r_kind - + real(r_double),parameter:: rmiss=10d7 ! Declare local variables logical outside,inflate_error - logical luse,ithinp + logical luse,ithinp,do_qc logical,allocatable,dimension(:,:):: lmsg ! set true when convinfo entry id found in a message character(70) obstr_v1, obstr_v2,hdrtr_v1,hdrtr_v2 @@ -172,7 +172,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis character(8) c_prvstg,c_sprvstg character(8) c_station_id,stationid - integer(i_kind) mxtb,nmsgmax + integer(i_kind) mxtb,nmsgmax,qcret integer(i_kind) ireadmg,ireadsb,iuse integer(i_kind) i,maxobs,idomsfc,nsattype,ncount integer(i_kind) nc,nx,isflg,itx,j,nchanl @@ -194,7 +194,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis integer(i_kind),dimension(5):: idate5 integer(i_kind),allocatable,dimension(:):: nrep,isort,iloc integer(i_kind),allocatable,dimension(:,:):: tab - + integer(i_kind) :: icnt(1000) integer(i_kind) ntime,itime @@ -265,6 +265,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis wjbmax=5.0_r_kind pflag=0 var_jb=zero + icnt=0 ! allocate(etabl(302,33,6)) ! add 2 ObsErr profiles for GOES-R IR(itype=301) and WV(itype=300) (not used yet, 2015-07-08, Genkova) @@ -366,6 +367,21 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis itype=254 endif endif + + else if(trim(subset) == 'NC005041' .or. trim(subset) == 'NC005042' .or. & + trim(subset) == 'NC005043') then + if( hdrdat(1) >=r100 .and. hdrdat(1) <=r199 ) then ! the range of JMA satellite IDS + if(hdrdat(9) == one) then ! IR winds + itype=252 + else if(hdrdat(9) == two) then ! visible winds + itype=242 + else if(hdrdat(9) == three) then ! WV cloud top + itype=250 + else if(hdrdat(9) >= four) then ! WV deep layer,monitored + itype=250 + endif + endif + else if(trim(subset) == 'NC005044' .or. trim(subset) == 'NC005045' .or. & trim(subset) == 'NC005046') then if( hdrdat(1) >=r100 .and. hdrdat(1) <=r199 ) then ! the range of JMA satellite IDS @@ -392,6 +408,25 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis itype=250 endif endif + + else if(trim(subset) == 'NC005001' .or. trim(subset) == 'NC005002' .or. & + trim(subset) == 'NC005003' ) then + if( hdrdat(1) >=r250 .and. hdrdat(1) <=r299 ) then ! the range of NESDIS satellite IDS + if(hdrdat(9) == one) then ! IR winds + if(hdrdat(12) <50000000000000.0_r_kind) then + itype=245 + else + itype=240 ! short wave IR winds + endif + else if(hdrdat(9) == two ) then ! visible winds + itype=251 + else if(hdrdat(9) == three ) then ! WV cloud top + itype=246 + else if(hdrdat(9) >= four ) then ! WV deep layer,monitored + itype=247 + endif + endif + else if(trim(subset) == 'NC005010' .or. trim(subset) == 'NC005011' .or. & trim(subset) == 'NC005012' ) then if( hdrdat(1) >=r250 .and. hdrdat(1) <=r299 ) then ! the range of NESDIS satellite IDS @@ -566,6 +601,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis !! read satellite winds one type a time loop_convinfo: do nx=1,ntread + ! set parameters for processing the next satwind type use_all = .true. use_all_tm = .true. ithin=0 @@ -607,8 +643,9 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis endif endif + ! Open and read the file once for each satwnd type call closbf(lunin) - close(lunin) + !close(lunin) open(lunin,file=trim(infile),form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) @@ -642,6 +679,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ee=r110 qifn=r110 qify=r110 + qm=2 ! Test for BUFR version using lat/lon mnemonics call ufbint(lunin,hdrdat_test,2,1,iret, 'CLAT CLON') @@ -653,18 +691,18 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis call ufbint(lunin,obsdat,4,1,iret,obstr_v1) endif + ! reject data with missing pressure or wind ppb=obsdat(2) - if (ppb > 100000000.0_r_kind .or. & - hdrdat(3) >100000000.0_r_kind .or. & - obsdat(4) > 100000000.0_r_kind) cycle loop_readsb - if(ppb >r10000) ppb=ppb/r100 + if(ppb>rmiss .or. hdrdat(3)>rmiss .or. obsdat(4)>rmiss) cycle loop_readsb + if(ppb>r10000) ppb=ppb/r100 ! ppb<10000 may indicate data reported in daPa or hPa + + ! reject date above 125mb (or 850 for regional) if (ppb r90 ) cycle loop_readsb if( hdrdat(3) =240.and.itype<=279) icnt(itype)=icnt(itype)+1 + + ! test for QCSTR or MANDATORY QC - if not skip over the extra blocks + call ufbrep(lunin,qcdat,3,12,qcret,qcstr) + do_qc = subset(1:7)=='NC00503'.and.nint(hdrdat(1))>=270 + do_qc = do_qc.or.subset(1:7)=='NC00501' + do_qc = do_qc.or.subset=='NC005081'.or.subset=='NC005091' + do_qc = do_qc.or.qcret>0 + + ! assign types and get quality info: start + + if(.not.do_qc) then + continue + else if(trim(subset) == 'NC005064' .or. trim(subset) == 'NC005065' .or. & trim(subset) == 'NC005066') then if( hdrdat(1) = r50) then ! the range of EUMETSAT satellite IDs c_prvstg='EUMETSAT' @@ -1000,7 +1052,7 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis endif ! get quality information THIS SECTION NEEDS TO BE TESTED!!! call ufbint(lunin,rep_array,1,1,iret, '{AMVIVR}') - irep_array = int(rep_array) + irep_array = max(1,int(rep_array)) allocate( amvivr(2,irep_array)) call ufbrep(lunin,amvivr,2,irep_array,iret, 'TCOV CVWD') pct1 = amvivr(2,1) ! use of pct1 (a new variable in the BUFR) is introduced by Nebuda/Genkova @@ -1259,17 +1311,17 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis if(itype==251 ) then; c_prvstg='GOESR' ; c_sprvstg='VIS' ; endif if(itype==241 ) then; c_prvstg='GOESR' ; c_sprvstg='IR' ; endif !to be revisited I.Genkova - endif + endif ! Extra block for GOES-R winds: End else ! wind is not recognised and itype is not assigned write(6,*) 'read_satwnd: WIND IS NOT RECOGNIZEd and we are in hell' cycle loop_readsb endif - if ( itype == -1 ) cycle loop_readsb ! unassigned itype - ! assign types and get quality info : end + if ( itype == -1 ) cycle loop_readsb ! unassigned itype + if ( qify == zero) qify=r110 if ( qifn == zero) qifn=r110 if ( ee == zero) ee=r110 diff --git a/src/gsi/read_ssmi.f90 b/src/gsi/read_ssmi.f90 index 3e47ee79b5..cece78ac03 100755 --- a/src/gsi/read_ssmi.f90 +++ b/src/gsi/read_ssmi.f90 @@ -132,7 +132,8 @@ subroutine read_ssmi(mype,val_ssmi,ithin,rmesh,jsatid,gstime,& real(r_kind),parameter:: tbmin=70.0_r_kind real(r_kind),parameter:: tbmax=320.0_r_kind character(80),parameter:: hdr1b='SAID YEAR MNTH DAYS HOUR MINU SECO ORBN' !use for ufbint() - character(40),parameter:: str1='CLAT CLON SFTG POSN SAZA' !use for ufbint() + character(40),parameter:: str1='CLATH CLONH SFTG POSN SAZA' !use for ufbint() new + character(40),parameter:: strx='CLAT CLON SFTG POSN SAZA' !use for ufbint() old character(40),parameter:: str2='TMBR' !use for ufbrep() ! Declare local variables @@ -302,6 +303,7 @@ subroutine read_ssmi(mype,val_ssmi,ithin,rmesh,jsatid,gstime,& ! SSM/I data are stored in groups of nscan, hence the loop. call ufbint(lnbufr,midat,nloc,nscan,iret,str1) + if(midat(1,1)>10e8) call ufbint(lnbufr,midat,nloc,nscan,iret,strx) !--- Extract brightness temperature data. Apply gross check to data. @@ -309,7 +311,6 @@ subroutine read_ssmi(mype,val_ssmi,ithin,rmesh,jsatid,gstime,& call ufbrep(lnbufr,mirad,1,nchanl*nscan,iret,str2) - ij=0 scan_loop: do js=1,nscan @@ -511,7 +512,6 @@ 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/setupoz.f90 b/src/gsi/setupoz.f90 index 24381df447..2008f37559 100644 --- a/src/gsi/setupoz.f90 +++ b/src/gsi/setupoz.f90 @@ -83,6 +83,9 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& ! 2017-02-09 guo - Remove m_alloc, n_alloc. ! . Remove my_node with corrected typecast(). ! 2017-10-27 todling - revised netcdf output for lay case; obs-sens needs attention +! 2020-02-26 todling - reset obsbin from hr to min +! 2022-08-10 karpowicz - fixes to ncdiag air_pressure_levels, change mass output to +! ppmv/mole fraction, fix ompsnm scan positoin and solar zenith angle. ! ! input argument list: ! lunin - unit from which to read observations @@ -110,10 +113,10 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& use mpeu_util, only: die,perr,getindex use kinds, only: r_kind,r_single,i_kind - use state_vectors, only: svars3d, levels + use state_vectors, only: svars3d, levels, nsdim use constants, only : zero,half,one,two,tiny_r_kind - use constants, only : rozcon,cg_term,wgtlim,h300,r10 + use constants, only : constoz,rozcon,cg_term,wgtlim,h300,r10,r100,r1000 use m_obsdiagNode, only : obs_diag use m_obsdiagNode, only : obs_diags @@ -131,6 +134,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& use m_obsLList, only : obsLList use obsmod, only : nloz_omi use obsmod, only : luse_obsdiag +! use obsmod, only : wrtgeovals use obsmod, only: netcdf_diag, binary_diag, dirname use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, & @@ -155,8 +159,8 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& implicit none ! !INPUT PARAMETERS: - type(obsLList ),target,dimension(:),intent(in):: obsLL - type(obs_diags),target,dimension(:),intent(in):: odiagLL + type(obsLList ),target,dimension(:),intent(inout):: obsLL + type(obs_diags),target,dimension(:),intent(inout):: odiagLL integer(i_kind) , intent(in ) :: lunin ! unit from which to read observations integer(i_kind) , intent(in ) :: mype ! mpi task id @@ -191,9 +195,9 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& ! Declare local variables - real(r_kind) omg,rat_err2,dlat,dtime,dlon,rat_err4diag + real(r_kind) omg,rat_err2,dlat,dtime,dlon real(r_kind) cg_oz,wgross,wnotgross,wgt,arg,exp_arg,term - real(r_kind) psi,errorinv + real(r_kind) psi,errorinv,rat_err4diag real(r_kind),dimension(nlevs):: ozges,varinv3,ozone_inv,ozobs,varinv4diag real(r_kind),dimension(nlevs):: ratio_errors,error real(r_kind),dimension(nlevs-1):: ozp @@ -201,6 +205,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& real(r_kind),dimension(nlevs):: pobs,gross,tnoise real(r_kind),dimension(nreal+nlevs,nobs):: data real(r_kind),dimension(nsig+1)::prsitmp + real(r_kind),dimension(nsig)::ozgestmp real(r_single),dimension(nlevs):: pob4,grs4,err4 real(r_single),dimension(ireal,nobs):: diagbuf real(r_single),allocatable,dimension(:,:,:)::rdiagbuf @@ -212,9 +217,10 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& real(r_kind),dimension(nsig,nloz_omi+1):: doz_dz1 integer(i_kind) :: oz_ind, nind, nnz type(sparr2) :: dhx_dx + real(r_single), dimension(nsdim) :: dhx_dx_array integer(i_kind) i,nlev,ii,jj,iextra,ibin, kk, nperobs - integer(i_kind) k,j,nz,jc,idia,irdim1,istatus,ioff0 + integer(i_kind) k1,k2,k,j,nz,jc,idia,irdim1,istatus,ioff0,ioff1 integer(i_kind) ioff,itoss,ikeep,ierror_toq,ierror_poq integer(i_kind) isolz,ifovn,itoqf integer(i_kind) mm1,itime,ilat,ilon,ilate,ilone,itoq,ipoq @@ -280,7 +286,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& iouse(jc)=iuse_oz(j) tnoise(jc)=error_oz(j) gross(jc)=min(r10*gross_oz(j),h300) - if (obstype == 'sbuv2' .or. obstype == 'ompsnp') then + if (obstype == 'sbuv2' .or. obstype == 'ompsnp' .or. obstype == 'ompsnpnc') then pobs(jc)=pob_oz(j) * 1.01325_r_kind else pobs(jc)=pob_oz(j) @@ -319,6 +325,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& irdim1=7 ioff0=irdim1 if(lobsdiagsave) irdim1=irdim1+4*miter+1 + ioff1=irdim1 if (save_jacobian) then nnz = nsig ! number of non-zero elements in dH(x)/dx profile nind = 1 @@ -363,7 +370,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& dlon=data(ilon,i) dtime=data(itime,i) - if (obstype == 'sbuv2' .or. obstype == 'ompsnp') then + if (obstype == 'sbuv2' .or. obstype == 'ompsnp' .or. obstype == 'ompsnpnc') then if (nobskeep>0) then ! write(6,*)'setupozlay: nobskeep',nobskeep call stop2(259) @@ -388,7 +395,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& enddo end if - if (obstype == 'omieff' .or. obstype == 'tomseff') then + if (obstype == 'omieff' .or. obstype == 'tomseff' .or. obstype == 'ompsnmeff') then pob_oz_omi(nloz_omi) = 1000.0_r_kind* 1.01325_r_kind do j=nloz_omi-1, 1, -1 pob_oz_omi(j) = pob_oz_omi(j+1)/2.0 @@ -410,7 +417,16 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& call grdcrd1(ozp_omi(nloz_omi),prsitmp,nsig+1,-1) end if - if (obstype /= 'omieff' .and. obstype /= 'tomseff') then + call tintrp2a1(ges_oz,ozgestmp,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + + ! need call to get pressures for pressure level output in ncdiags + call tintrp2a1(ges_prsi,prsitmp,dlat,dlon,dtime,hrdifsig,& + nsig+1,mype,nfldsig) + + + if (obstype /= 'omieff' .and. obstype /= 'tomseff' .and. & + obstype /= 'ompsnmeff' ) then call intrp3oz1(ges_oz,ozges,dlat,dlon,ozp,dtime,& nlevs,mype,doz_dz) endif @@ -441,7 +457,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& ! For OMI/GOME, nlev=1 do k=1,nlev j=ipos(k) - if (obstype == 'omieff' .or. obstype == 'tomseff' ) then + if (obstype == 'omieff' .or. obstype == 'tomseff' .or. obstype == 'ompsnmeff') then ioff=ifovn+1 ! else ioff=nreal+k ! SBUV and OMI w/o efficiency factors @@ -449,7 +465,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& ! Compute innovation and load obs error into local array ! KW OMI and TOMS have averaging kernels - if (obstype == 'omieff' .or. obstype == 'tomseff' ) then + if (obstype == 'omieff' .or. obstype == 'tomseff' .or. obstype == 'ompsnmeff') then ! everything in data is from top to bottom nlayers = nloz_omi + 1 apriori(1:nloz_omi) = data(ioff:ioff+nloz_omi -1, i) @@ -542,7 +558,8 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& rdiagbuf(3,k,ii) = errorinv ! inverse observation error if (obstype == 'gome' .or. obstype == 'omieff' .or. & obstype == 'omi' .or. obstype == 'tomseff' .or. & - obstype == 'ompstc8') then + obstype == 'ompsnmeff' .or. obstype == 'ompstc8' .or. & + obstype == 'ompsnm') then rdiagbuf(4,k,ii) = data(isolz,i) ! solar zenith angle rdiagbuf(5,k,ii) = data(ifovn,i) ! field of view number else @@ -556,7 +573,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& endif rdiagbuf(7,k,ii) = 1.e+10_r_single ! spread (filled in by EnKF) - idia = ioff0 + idia = ioff1 if (save_jacobian) then oz_ind = getindex(svars3d, 'oz') if (oz_ind < 0) then @@ -574,18 +591,39 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& endif if (netcdf_diag) then + k1 = k + k2 = k - 1 + if(k2 == 0)k2 = 1 + if(k == nlevs)then + k1=nlevs-1 + k2=1 + endif + if (obstype == 'sbuv2' .or. obstype == 'ompsnp' .or. obstype == 'ompsnpnc' ) then + call nc_diag_metadata("TopLevelPressure",sngl(pobs(k2)*r100)) + call nc_diag_metadata("BottomLevelPressure", & + sngl(pobs(k1)*r100)) + else + call & + nc_diag_metadata("TopLevelPressure",sngl(prsitmp(nsig+1)*r1000) ) + call nc_diag_metadata("BottomLevelPressure", & + sngl(prsitmp(1)*r1000) ) + endif call nc_diag_metadata("MPI_Task_Number", mype ) call nc_diag_metadata("Latitude", sngl(data(ilate,i)) ) call nc_diag_metadata("Longitude", sngl(data(ilone,i)) ) call nc_diag_metadata("Time", sngl(data(itime,i)-time_offset) ) - call nc_diag_metadata("Reference_Pressure", sngl(pobs(k)) ) + call nc_diag_metadata("Reference_Pressure",sngl(pobs(k)*r100) ) call nc_diag_metadata("Analysis_Use_Flag", iouse(k) ) call nc_diag_metadata("Observation", sngl(ozobs(k))) call nc_diag_metadata("Inverse_Observation_Error", sngl(errorinv)) + call nc_diag_metadata("Input_Observation_Error", sngl(error(k))) call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(ozone_inv(k))) call nc_diag_metadata("Obs_Minus_Forecast_unadjusted",sngl(ozone_inv(k))) + call nc_diag_metadata("Forecast_unadjusted", sngl(ozges(k))) + call nc_diag_metadata("Forecast_adjusted",sngl(ozges(k))) if (obstype == 'gome' .or. obstype == 'omieff' .or. & - obstype == 'omi' .or. obstype == 'tomseff' ) then + obstype == 'omi' .or. obstype == 'tomseff' .or. & + obstype == 'ompsnmeff' .or. obstype == 'ompsnm') then call nc_diag_metadata("Solar_Zenith_Angle", sngl(data(isolz,i)) ) call nc_diag_metadata("Scan_Position", sngl(data(ifovn,i)) ) else @@ -598,10 +636,13 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& call nc_diag_metadata("Row_Anomaly_Index", sngl(rmiss) ) endif if (save_jacobian) then - call nc_diag_data2d("Observation_Operator_Jacobian_stind", dhx_dx%st_ind) - call nc_diag_data2d("Observation_Operator_Jacobian_endind", dhx_dx%end_ind) - call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single)) + call fullarray(dhx_dx, dhx_dx_array) + call nc_diag_data2d("Observation_Operator_Jacobian", dhx_dx_array) endif + !if (wrtgeovals) then + ! call nc_diag_data2d("mole_fraction_of_ozone_in_air", sngl(constoz*ozgestmp)) + ! call nc_diag_data2d("air_pressure_levels",sngl(prsitmp*r1000)) + !endif endif endif @@ -642,7 +683,8 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& my_head%elon= data(ilone,i) nlevp=max(nlev-1,1) - if (obstype == 'omieff' .or. obstype == 'tomseff' ) nlevp = nloz_omi + if (obstype == 'omieff' .or. obstype == 'tomseff' .or. & + obstype == 'ompsnmeff') nlevp = nloz_omi allocate(my_head%res(nlev), & my_head%err2(nlev), & my_head%raterr2(nlev), & @@ -679,11 +721,12 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& my_head%luse=luse(i) my_head%time=dtime - if (obstype == 'sbuv2'.or. obstype == 'ompsnp' ) then + if (obstype == 'sbuv2' .or. obstype == 'ompsnp' .or. obstype == 'ompsnpnc') then do k=1,nlevs-1 my_head%prs(k) = ozp(k) enddo - else if (obstype == 'omieff' .or. obstype == 'tomseff') then + else if (obstype == 'omieff' .or. obstype == 'tomseff' .or. & + obstype == 'ompsnmeff') then do k=1,nloz_omi my_head%prs(k) = ozp_omi(k) enddo @@ -742,7 +785,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& if (ozone_diagsave.and.lobsdiagsave.and.luse(i)) then associate(odiag => my_diag) - idia=6 + idia=ioff0 do jj=1,miter idia=idia+1 if (odiag%muse(jj)) then @@ -785,7 +828,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& if(in_curbin) then if (ozone_diagsave.and.lobsdiagsave.and.luse(i)) then - rdiagbuf(7:irdim1,1:nlevs,ii) = zero + rdiagbuf(ioff0+1:irdim1,1:nlevs,ii) = zero endif endif ! (in_curbin) @@ -932,6 +975,7 @@ subroutine init_netcdf_diag_ call nc_diag_header("Satellite_Sensor", isis) call nc_diag_header("Satellite", dplat(is)) call nc_diag_header("Observation_type", obstype) + call nc_diag_header("Number_of_state_vars", nsdim ) call nc_diag_header("pobs", pobs) call nc_diag_header("gross",gross) call nc_diag_header("tnoise",tnoise) @@ -995,7 +1039,9 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& ! 2016-12-09 mccarty - add netcdf_diag capability ! 2017-02-09 guo - Remove m_alloc, n_alloc. ! . Remove my_node with corrected typecast(). -! +! 2020-02-26 todling - reset obsbin from hr to min +! 2022-08-10 karpowicz/todling - replace ncdiag analysis use flag with +/-1 instead of zero +! ! input argument list: ! lunin - unit from which to read observations ! mype - mpi task id @@ -1036,6 +1082,7 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& use obsmod, only : mype_diaghdr,dirname,time_offset,ianldate use obsmod, only : lobsdiag_allocated,lobsdiagsave,lobsdiag_forenkf use obsmod, only: netcdf_diag, binary_diag, dirname +! use obsmod, only: wrtgeovals use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, & nc_diag_write, nc_diag_data2d use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close @@ -1048,7 +1095,7 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& use guess_grids, only : nfldsig,ges_lnprsl,hrdifsig use constants, only : zero,half,one,two,tiny_r_kind,four - use constants, only : cg_term,wgtlim,r10,constoz + use constants, only : cg_term,wgtlim,r10,r100,r1000,constoz use gsi_4dvar, only: nobs_bins,hr_obsbin @@ -1066,8 +1113,8 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& implicit none ! !INPUT PARAMETERS: - type(obsLList ),target,dimension(:),intent(in):: obsLL - type(obs_diags),target,dimension(:),intent(in):: odiagLL + type(obsLList ),target,dimension(:),intent(inout):: obsLL + type(obs_diags),target,dimension(:),intent(inout):: odiagLL integer(i_kind) , intent(in ) :: lunin ! unit from which to read observations integer(i_kind) , intent(in ) :: mype ! mpi task id @@ -1106,26 +1153,30 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& real(r_kind) o3ges, o3ppmv real(r_kind) rlow,rhgh,sfcchk - real(r_kind) omg,rat_err2,dlat,dtime,dlon,rat_err4diag + real(r_kind) omg,rat_err2,dlat,dtime,dlon real(r_kind) cg_oz,wgross,wnotgross,wgt,arg,exp_arg,term real(r_kind) errorinv - real(r_kind) psges,ozlv,airnd,uvnd,visnd + real(r_kind) psges,ozlv, airnd, uvnd, visnd - real(r_kind) varinv3,ratio_errors,varinv4diag + real(r_kind) varinv3,ratio_errors + real(r_kind) varinv4diag,rat_err4diag real(r_kind) dpres,obserror,ozone_inv,preso3l real(r_kind),dimension(nreal+nlevs,nobs):: data real(r_kind),dimension(nsig):: prsltmp real(r_single),dimension(ireal,nobs):: diagbuf real(r_single),allocatable,dimension(:,:,:)::rdiagbuf + real(r_kind),dimension(nsig+1)::prsitmp + real(r_kind),dimension(nsig)::ozgestmp integer(i_kind) i,ii,jj,iextra,ibin - integer(i_kind) k,j,idia,irdim1,ioff0 + integer(i_kind) k1,k2,k,j,idia,irdim1,ioff0,ioff1 integer(i_kind) isolz,iuse integer(i_kind) mm1,itime,ilat,ilon,ilate,ilone,iozmr,ilev,ipres,iprcs,imls_levs + integer(i_kind) iairnd, iuvnd, ivisnd integer(i_kind),dimension(iint,nobs):: idiagbuf - integer(i_kind) iairnd,iuvnd,ivisnd real(r_kind) gross,tnoise,pobs + character(12) string character(10) filex character(128) diag_ozone_file @@ -1178,6 +1229,7 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& irdim1=10 ioff0 = irdim1 if(lobsdiagsave) irdim1=irdim1+4*miter+1 + ioff1 = irdim1 if (save_jacobian) then nnz = 2 ! number of non-zero elements in dH(x)/dx profile nind = 1 @@ -1315,6 +1367,9 @@ subroutine setupozlev(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& obserror=1.0e6_r_kind end if + call tintrp2a1(ges_oz,ozgestmp,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + ! Interpolate guess ozone to observation location and time call tintrp31(ges_oz,o3ges,dlat,dlon,dpres,dtime, & hrdifsig,mype,nfldsig) @@ -1627,7 +1682,7 @@ subroutine contents_binary_diag_(odiag) rdiagbuf(10,1,ii) = visnd ! log10 ozone number density vis if (lobsdiagsave) then - idia=6 + idia=ioff0 do jj=1,miter idia=idia+1 if (odiag%muse(jj)) then @@ -1649,6 +1704,7 @@ subroutine contents_binary_diag_(odiag) rdiagbuf(idia,1,ii) = odiag%obssen(jj) enddo endif + idia = ioff1 if (save_jacobian) then call writearray(dhx_dx, rdiagbuf(idia+1:irdim1,1,ii)) idia = idia + size(dhx_dx) @@ -1668,19 +1724,33 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Observation", sngl(ozlv) ) call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(ozone_inv) ) call nc_diag_metadata("Obs_Minus_Forecast_unadjusted",sngl(ozone_inv) ) - call nc_diag_metadata("Reference_Pressure", sngl(preso3l) ) + call nc_diag_metadata("Reference_Pressure", sngl(preso3l*r100) ) ! Pa + if(luse(i)) then + call nc_diag_metadata("Analysis_Use_Flag", one ) + else + call nc_diag_metadata("Analysis_Use_Flag", -one ) + endif call nc_diag_metadata("Input_Observation_Error", sngl(obserror) ) - if(obstype =="omps_lp")then + if(obstype =="ompslp")then call nc_diag_metadata("Log10 Air Number Density", sngl(airnd)) - call nc_diag_metadata("Log10 Ozone Number Density UV", sngl(uvnd)) + call nc_diag_metadata("Log10 Ozone Number Density UV", sngl(uvnd)) call nc_diag_metadata("Log10 Ozone Number Density VIS",sngl(visnd)) endif - - if(luse(i)) then - call nc_diag_metadata("Analysis_Use_Flag", 1 ) - else - call nc_diag_metadata("Analysis_Use_Flag", -1 ) + call nc_diag_metadata("Forecast_adjusted", sngl(o3ppmv)) + call nc_diag_metadata("Forecast_unadjusted", sngl(o3ppmv)) + !if (wrtgeovals) then + ! ozgestmp = ozgestmp *constoz + ! call nc_diag_data2d("mole_fraction_of_ozone_in_air", sngl(ozgestmp)) + ! call nc_diag_data2d("air_pressure",sngl(exp(prsltmp)*r1000)) ! Pa + !endif + k1 = k + k2 = k - 1 + if(k2 == 0)k2 = 1 + if(k == nlevs)then + k1=nlevs-1 + k2=1 endif + if (save_jacobian) then call nc_diag_data2d("Observation_Operator_Jacobian_stind", dhx_dx%st_ind) call nc_diag_data2d("Observation_Operator_Jacobian_endind", dhx_dx%end_ind)