Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 34 additions & 15 deletions src/gsi/read_prepbufr.F90
Original file line number Diff line number Diff line change
Expand Up @@ -340,7 +340,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
real(r_kind) qtflg,tdry,rmesh,ediff,usage,ediff_ps,ediff_q,ediff_t,ediff_uv,ediff_pw
real(r_kind) u0,v0,uob,vob,rgustob,dx,dy,dx1,dy1,w00,w10,w01,w11
real(r_kind) qoe,qobcon,pwoe,pwmerr,dlnpob,ppb,poe,gustoe,visoe,qmaxerr
real(r_kind) toe,woe,errout,oelev,dlat,dlon,sstoe,dlat_earth,dlon_earth
real(r_kind) toe,woe,errout,oelev,dlat,dlon,sstoe,dlat_earth,dlon_earth, t29
real(r_kind) tdoe,mxtmoe,mitmoe,pmoe,howvoe,cldchoe
real(r_kind) dlat_earth_deg,dlon_earth_deg
real(r_kind) selev,elev,stnelev
Expand Down Expand Up @@ -519,7 +519,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
iqm = 11
iuse = 13
else if(psob) then
nreal=20
nreal=21
iqm=10
iuse = 12
else if(qob) then
Expand Down Expand Up @@ -724,6 +724,8 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&

! Extract type information
call ufbint(lunin,hdr,4,1,iret,hdstr2)
t29 = hdr(3)
it29 = nint(t29,r_double)
Comment thread
RussTreadon-NOAA marked this conversation as resolved.
kx=hdr(1)
if (aircraft_t_bc .and. acft_profl_file) then
kx0=kx
Expand Down Expand Up @@ -766,7 +768,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&

! identify drifting buoys - TYP=180/280 T29=562 and last three digits of SID between 500 and 999
! (see https://www.wmo.int/pages/prog/amp/mmop/wmo-number-rules.html) Set kx to 199/299
if (id_drifter .and. (kx==180 .or. kx==280) .and. nint(hdr(3),r_double)==562) then
if (id_drifter .and. (kx==180 .or. kx==280) .and. it29==562) then
rstation_id=hdr(4)
read(c_station_id,*,iostat=ios) iwmo
if (ios == 0 .and. iwmo > 0) then
Expand All @@ -776,7 +778,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
end if
end if

if (id_ship .and. (kx==180) .and. (nint(hdr(3),r_double)==522 .or. nint(hdr(3),r_double)==523)) then
if (id_ship .and. (kx==180) .and. (it29==522 .or. it29==523)) then
rstation_id=hdr(4)
kx = kx + 18
end if
Expand All @@ -793,15 +795,14 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&

! temporary specify iobsub until put in bufr file
iobsub = 0
if(kx == 280 .or. kx == 180 ) iobsub=hdr(3)
if(kx == 280 .or. kx ==180) then
if ( hdr(3) >555.0_r_kind .and. hdr(3) <565.0_r_kind ) then
if ( it29 > 555 .and. it29 < 565.0_r_kind ) then
iobsub=00
else
iobsub=01
endif
! Set saildrone to subtype 02
if (nint(hdr(3)) == 560) iobsub = 02
if (it29 == 560) iobsub = 02
endif
! Su suggested to keep both 289 and 290. But trunk only keep 290
! if(kx == 289 .or. kx == 290) iobsub=hdr(2)
Expand Down Expand Up @@ -1036,6 +1037,8 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
! Extract type, date, and location information
call ufbint(lunin,hdr,8,1,iret,hdstr)
kx=hdr(5)
t29 = hdr(8)
it29 = nint(t29,r_double)

if (.not.(aircraft_t_bc .and. acft_profl_file)) then
if(abs(hdr(3))>r90 .or. abs(hdr(2))>r360) cycle loop_readsb
Expand All @@ -1050,7 +1053,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
!
! identify drifting buoys - TYP=180/280 T29=562 and last three digits of SID between 500 and 999
! (see https://www.wmo.int/pages/prog/amp/mmop/wmo-number-rules.html) Set kx to 199/299
if (id_drifter .and. (kx==180 .or. kx==280) .and. nint(hdr(8),r_double)==562) then
if (id_drifter .and. (kx==180 .or. kx==280) .and. it29==562) then
rstation_id=hdr(1)
read(c_station_id,*,iostat=ios) iwmo
if (ios == 0 .and. iwmo > 0) then
Expand All @@ -1060,7 +1063,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
end if
end if

if (id_ship .and. (kx==180) .and. (nint(hdr(8),r_double)==522 .or. nint(hdr(8),r_double)==523) ) then
if (id_ship .and. (kx==180) .and. (it29==522 .or. it29==523) ) then
rstation_id=hdr(1)
kx = kx + 18
end if
Expand Down Expand Up @@ -2147,7 +2150,6 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
oelev=r10+selev
endif
if (kx == 280 )then
it29=nint(hdr(8))
if(it29 == 522 .or. it29 == 523 .or. it29 == 531)then
! oelev=r20+selev
oelev=r20
Expand Down Expand Up @@ -2312,7 +2314,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
cdata_all(3,iout)=dlat ! grid relative latitude
cdata_all(4,iout)=dlnpob ! ln(pressure in cb)

cdata_all(5,iout)=obsdat(3,k)+t0c ! temperature ob.
cdata_all(5,iout)=obsdat(3,k)+t0c ! temperature ob.
cdata_all(6,iout)=rstation_id ! station id
cdata_all(7,iout)=t4dv ! time
cdata_all(8,iout)=nc ! type
Expand All @@ -2332,7 +2334,11 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
cdata_all(21,iout)=zz ! terrain height at ob location
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(24,iout)=obsdat(10,k) ! cat
if (it29 > 500) then
cdata_all(24,iout)=t29 ! it29
else
cdata_all(24,iout)=obsdat(10,k) ! cat
end if
cdata_all(25,iout)=var_jb(3,k) ! non linear qc for T
if (aircraft_t_bc_pof .or. aircraft_t_bc .or.aircraft_t_bc_ext) then
cdata_all(26,iout)=aircraftwk(1,k) ! phase of flight
Expand Down Expand Up @@ -2398,7 +2404,11 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
cdata_all(21,iout)=zz ! terrain height at ob location
cdata_all(22,iout)=r_prvstg(1,1) ! provider name
cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(24,iout)=obsdat(10,k) ! cat
if (it29 > 500) then
cdata_all(24,iout)=t29 ! it29
else
cdata_all(24,iout)=obsdat(10,k) ! cat
end if
cdata_all(25,iout)=var_jb(5,k) ! non linear qc parameter
cdata_all(26,iout)=one ! hilbert curve weight, modified later
cdata_all(27,iout)=windbiasfact ! bias correction factor
Expand Down Expand Up @@ -2472,7 +2482,12 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
cdata_all(18,iout)=r_prvstg(1,1) ! provider name
cdata_all(19,iout)=r_sprvstg(1,1) ! 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 (it29 > 500) then
cdata_all(21,iout)=t29 ! it29
else
cdata_all(21,iout)=obsdat(10,k) ! cat
end if
if(perturb_obs)cdata_all(22,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))

Expand Down Expand Up @@ -2521,7 +2536,11 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
cdata_all(19,iout)=zz ! terrain height at ob location
cdata_all(20,iout)=r_prvstg(1,1) ! provider name
cdata_all(21,iout)=r_sprvstg(1,1) ! subprovider name
cdata_all(22,iout)=obsdat(10,k) ! cat
if (it29 > 500) then
cdata_all(22,iout)=t29 ! it29
else
cdata_all(22,iout)=obsdat(10,k) ! cat
end if
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) &
Expand Down
6 changes: 4 additions & 2 deletions src/gsi/setupps.f90
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
real(r_kind),dimension(nele,nobs):: data
real(r_single),allocatable,dimension(:,:)::rdiagbuf

integer(i_kind) ier,ilon,ilat,ipres,ihgt,itemp,id,itime,ikx,iqc,iptrb,ijb
integer(i_kind) ier,ilon,ilat,ipres,ihgt,itemp,id,itime,ikx,iqc,iptrb,ijb,icat
integer(i_kind) ier2,iuse,ilate,ilone,istnelv,idomsfc,izz,iprvd,isprvd
integer(i_kind) ikxx,nn,ibin,ioff,ioff0
integer(i_kind) i,j,nchar,nreal,ii,jj,k,l,mm1
Expand Down Expand Up @@ -249,7 +249,8 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
iprvd=18 ! index of observation provider
isprvd=19 ! index of observation subprovider
ijb=20 ! index of non linear qc parameter
iptrb=21 ! index of ps perturbation
icat=21 ! index of ps perturbation
iptrb=22 ! index of ps perturbation

! Declare local constants
halfpi = half*pi
Expand Down Expand Up @@ -890,6 +891,7 @@ subroutine contents_netcdf_diag_(odiag)
call nc_diag_metadata("Observation_Class", obsclass )
call nc_diag_metadata("Observation_Type", ictype(ikx) )
call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) )
call nc_diag_metadata_to_single("Observation_Category", data(icat,i) )
!Replace direct calls to nc_diag_metadata with the screening subroutine
call nc_diag_metadata_to_single("Latitude", data(ilate,i) )
call nc_diag_metadata_to_single("Longitude", data(ilone,i) )
Expand Down
2 changes: 2 additions & 0 deletions src/gsi/setupq.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1380,6 +1380,7 @@ subroutine contents_netcdf_diag_(odiag)
call nc_diag_metadata("Observation_Class", obsclass )
call nc_diag_metadata("Observation_Type", ictype(ikx) )
call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) )
call nc_diag_metadata_to_single("Observation_Category", data(icat,i) )
call nc_diag_metadata_to_single("Latitude", data(ilate,i) )
call nc_diag_metadata_to_single("Longitude", data(ilone,i) )
! this is the obs height after being interpolated to the model (=model height)
Expand Down Expand Up @@ -1450,6 +1451,7 @@ subroutine contents_netcdf_diagp_(odiag)
call nc_diag_metadata("Observation_Class", obsclass )
call nc_diag_metadata("Observation_Type", ictype(ikx) )
call nc_diag_metadata("Observation_Subtype", -1 ) ! (-1 for pseudo obs sub-type)
call nc_diag_metadata_to_single("Observation_Category", data(icat,i) )
call nc_diag_metadata_to_single("Latitude", data(ilate,i) )
call nc_diag_metadata_to_single("Longitude", data(ilone,i) )
call nc_diag_metadata_to_single("Station_Elevation",data(istnelv,i) )
Expand Down
2 changes: 2 additions & 0 deletions src/gsi/setupt.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1823,6 +1823,7 @@ subroutine contents_netcdf_diag_(odiag)
call nc_diag_metadata("Observation_Class", obsclass )
call nc_diag_metadata("Observation_Type", ictype(ikx) )
call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) )
call nc_diag_metadata_to_single("Observation_Category", data(icat,i) )
call nc_diag_metadata_to_single("Latitude",data(ilate,i))
call nc_diag_metadata_to_single("Longitude",data(ilone,i))
! this is the obs height after being interpolated to the model (=model height)
Expand Down Expand Up @@ -1927,6 +1928,7 @@ subroutine contents_netcdf_diagp_(odiag)
call nc_diag_metadata("Observation_Class", obsclass )
call nc_diag_metadata("Observation_Type", ictype(ikx) )
call nc_diag_metadata("Observation_Subtype", -1 ) ! (-1 for pseudo obs sub-type)
call nc_diag_metadata_to_single("Observation_Category", data(icat,i) )
call nc_diag_metadata_to_single("Latitude",data(ilate,i))
call nc_diag_metadata_to_single("Longitude",data(ilone,i))
call nc_diag_metadata_to_single("Station_Elevation",data(istnelv,i))
Expand Down
1 change: 1 addition & 0 deletions src/gsi/setupw.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1782,6 +1782,7 @@ subroutine contents_netcdf_diag_(udiag,vdiag)
call nc_diag_metadata("Observation_Class", obsclass )
call nc_diag_metadata("Observation_Type", ictype(ikx) )
call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) )
call nc_diag_metadata_to_single("Observation_Category", data(icat,i) )
call nc_diag_metadata_to_single("Latitude",data(ilate,i) )
call nc_diag_metadata_to_single("Longitude",data(ilone,i) )
call nc_diag_metadata_to_single("Station_Elevation",data(ielev,i) )
Expand Down
14 changes: 7 additions & 7 deletions ush/prune_4nco_global.sh
Original file line number Diff line number Diff line change
Expand Up @@ -115,11 +115,11 @@ done

# Process ush directories and files
cd $topdir/ush
rlist="sub"
rlist="sub* run_observer"
for type in $rlist; do
if [[ "$mode" = "prune" ]]; then
if [ -e $type ]; then
git $string ${type}*
if [ -e ${type} ]; then
git $string ${type}
rc=$?
if [[ $rc -ne 0 ]]; then
echo "***ERROR*** git $string ${type}"
Expand All @@ -128,12 +128,12 @@ for type in $rlist; do
fi
elif [[ "$mode" = "restore" ]]; then
if [[ "$use_checkout" = "YES" ]]; then
git reset HEAD ${type}*
git checkout ${type}*
git reset HEAD ${type}
git checkout ${type}
rc=$?
else
git restore --staged ${type}*
git restore ${type}*
git restore --staged ${type}
git restore ${type}
rc=$?
fi
if [[ $rc -ne 0 ]]; then
Expand Down
27 changes: 27 additions & 0 deletions ush/run_observer/README.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
Script to run the GSI observer in standalone mode.
--------------------------------------------------

To get inital files from operations:

> sbatch rungsi_v17_dev.sh

You will need to set bdate in this script to the correct date.



To run the standalone code:

> sbatch rungsi_v17_dev.sh


This is basically a modified exglobal_atmos_analysis.sh.
To get started you will need to modify some of the following lines:
-----
# Path to previously built GSI is here
GSIDIR=/scratch3/NCEPDEV/da/${USER}/git/GSI_develop <<<<<<<<<<
machine=ursa <<<<<<<<<<

# Set experiment name and analysis date
adate=2025041500
exp=v17-dev.$adate

84 changes: 84 additions & 0 deletions ush/run_observer/get_initial_files_ops.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
#!/bin/sh
#SBATCH -o getfiles.out
#SBATCH -e getfiles.err
#SBATCH -p u1-service
#SBATCH --ntasks=1
#SBATCH --time=6:00:00
#SBATCH -A da-cpu
#SBATCH -J get_hpss_stuff

bdate=2025041500 # For stand alone GSI just need to change this date
edate=$bdate
cdate=${bdate}

expid=GDAS-ops
prod=v16.3

hpss_base_dir=/NCEPPROD/hpssprod/runhistory

NDATE=/home/Andrew.Collard/bin/ndate
Comment thread
RussTreadon-NOAA marked this conversation as resolved.
HPSSTAR=/home/Andrew.Collard/bin/hpsstar

while [[ ${cdate} -le ${edate} ]]; do

# for current analysis cycle
data_dir=/scratch3/NCEPDEV/stmp/$USER/$expid
mkdir -p ${data_dir}
cd ${data_dir}

y4a=`echo $cdate | cut -c1-4`
m2a=`echo $cdate | cut -c5-6`
d2a=`echo $cdate | cut -c7-8`
h2a=`echo $cdate | cut -c9-10`

yyyymmdda=${y4a}${m2a}${d2a}
hha=${h2a}
yyyymmddhha=${yyyymmdda}${hha}

# get required initial files from previous cycle
gdate=`$NDATE -6 ${cdate}`
y4g=`echo $gdate | cut -c1-4`
m2g=`echo $gdate | cut -c5-6`
d2g=`echo $gdate | cut -c7-8`
h2g=`echo $gdate | cut -c9-10`

yyyymmddg=${y4g}${m2g}${d2g}
hhg=${h2g}
yyyymmddhhg=${yyyymmddg}${hhg}

hpssa_dir=${hpss_base_dir}/rh${y4a}/${y4a}${m2a}/${yyyymmdda}
hpssg_dir=${hpss_base_dir}/rh${y4g}/${y4g}${m2g}/${yyyymmddg}

# =====================
# Get deterministics
# =====================
# Get initial conditions related to bias correction from previous (guess) cycle
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_restart.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.abias
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_restart.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.abias_pc
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_restart.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.abias_air
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_restart.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.radstat

# Get initial conditions related to model states from previous(guess) cycle
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_nc.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.atmf003.nc
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_nc.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.atmf006.nc
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_nc.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.atmf009.nc
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_nc.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.sfcf003.nc
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_nc.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.sfcf006.nc
$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_gdas.${yyyymmddg}_${h2g}.gdas_nc.tar ./gdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.sfcf009.nc

$HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_enkfgdas.${yyyymmddg}_${h2g}.enkfgdas.tar ./enkfgdas.${yyyymmddg}/${hhg}/atmos/gdas.t${hhg}z.sfcf006.ensmean.nc

# =============================================================
# Get ensembles
# ./enkfgdas.20200914/18/atmos/mem002/gdas.t18z.atmf003.nc
# =============================================================
# $HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_enkfgdas.${yyyymmddg}_${h2g}.enkfgdas_restart_grp1.tar . &
# $HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_enkfgdas.${yyyymmddg}_${h2g}.enkfgdas_restart_grp1.tar . &
# $HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_enkfgdas.${yyyymmddg}_${h2g}.enkfgdas_restart_grp1.tar . &
# $HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_enkfgdas.${yyyymmddg}_${h2g}.enkfgdas_restart_grp1.tar . &

# $HPSSTAR get ${hpssg_dir}/com_gfs_${prod}_enkfgdas.${yyyymmddg}_${h2g}.enkfgdas_restart_grp1.tar . &

cdate=`$NDATE +6 ${cdate}`
Comment thread
RussTreadon-NOAA marked this conversation as resolved.
done

Loading