From d199ed59e9dccdcafb1f2b44125c183376b0cdf9 Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Thu, 27 May 2021 19:08:39 +0000 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#120. merge GFS v16.1.1 DA changes into master --- doc/Release_Notes.gfsda.v16.1.0.txt | 121 ++++ fix | 2 +- scripts/exgdas_atmos_chgres_forenkf.sh | 10 +- scripts/exgdas_enkf_ecen.sh | 46 +- scripts/exgdas_enkf_fcst.sh | 12 +- scripts/exgdas_enkf_post.sh | 6 +- scripts/exgdas_enkf_select_obs.sh | 6 +- scripts/exgdas_enkf_sfc.sh | 11 +- scripts/exgdas_enkf_update.sh | 11 +- scripts/exglobal_atmos_analysis.sh | 34 +- scripts/exglobal_atmos_analysis_calc.sh | 12 +- scripts/exglobal_diag.sh | 20 +- src/enkf/gridio_gfs.f90 | 73 +- src/fv3gfs_ncio/module_fv3gfs_ncio.f90 | 2 +- src/gsi/cplr_gfs_ensmod.f90 | 8 +- src/gsi/gesinfo.F90 | 12 +- src/gsi/gsimod.F90 | 6 +- src/gsi/guess_grids.F90 | 5 +- src/gsi/ncepnems_io.f90 | 18 +- src/gsi/netcdfgfs_io.f90 | 26 +- src/gsi/read_files.f90 | 45 +- src/gsi/read_fl_hdob.f90 | 9 - src/gsi/setupbend.f90 | 671 +++++++++--------- src/gsi/setuprad.f90 | 2 +- src/gsi/stop1.f90 | 13 +- src/gsi/write_incr.f90 | 6 +- ush/build_all_cmake.sh | 23 +- ush/prune_4nco_global.sh | 189 +++++ .../scripts/exgdas_atmos_vminmon.sh | 1 - .../gfs.v1.0.0/scripts/exgfs_atmos_vminmon.sh | 1 - .../scripts/exgdas_atmos_verfrad.sh | 10 +- 31 files changed, 870 insertions(+), 541 deletions(-) create mode 100644 doc/Release_Notes.gfsda.v16.1.0.txt create mode 100755 ush/prune_4nco_global.sh diff --git a/doc/Release_Notes.gfsda.v16.1.0.txt b/doc/Release_Notes.gfsda.v16.1.0.txt new file mode 100644 index 0000000000..ca805b08ee --- /dev/null +++ b/doc/Release_Notes.gfsda.v16.1.0.txt @@ -0,0 +1,121 @@ +GFSDA.v16.1.0 RELEASE NOTES + + +PRELUDE + +* The tag gfsda.v16.1.0 was created in the NOAA-EMC/GSI repository to + support the GFSv16.1 implementation and is included as a part of the + NOAA-EMC/global-workflow tag EMC-v16.1.0. See release notes + docs/Release_Notes.gfs.v16.1.0.txt in the NOAA-EMC/global-workflow + tag for more details on the motivation, workflow changes, and + implementation instructions. + + +FIX CHANGES + +* fix/ + * global_convinfo.txt: Change GeoOptics (SAID=265) from -1 to 1 + for active assimilation + + +SCRIPTS CHANGES + +* scripts/ + * exgdas_atmos_chgres_forenkf.sh: Addresses GFSv16 bugzilla #1198 + * exgdas_enkf_ecen.sh: Addresses GFSv16 bugzilla #1198 + * exgdas_enkf_fcst.sh: Addresses GFSv16 bugzilla #1198 + * exgdas_enkf_post.sh: Addresses GFSv16 bugzilla #1198 + * exgdas_enkf_select_obs.sh: Addresses GFSv16 bugzilla #1198 + * exgdas_enkf_sfc.sh: Addresses GFSv16 bugzilla #1198 + * exgdas_enkf_update.sh: Addresses GFSv16 bugzilla #1198 + * exglobal_atmos_analysis.sh: Addresses GFSv16 bugzilla #1198, + adds default values for commercial RO parameters + * exglobal_atmos_analysis_calc.sh: Addresses GFSv16 bugzilla #1198 + * exglobal_diag.sh: Addresses GFSv16 bugzilla #1198 + + +SRC CHANGES + +* src/ + * enkf/ + * gridio_gfs.f90: Changes delz calculation to negative definite + + * fv3gfs_ncio/ + * module_fv3gfs_ncio.f90: Addresses GFSv16 bugzilla #1196 + + * gsi/ + * cplr_gfs_ensmod.f90: Addresses GFSv16 bugzilla #1196 + * genstats_gps.f90: GPSRO bugfix + * gesinfo.F90: Addresses GFSv16 bugzilla #1196 + * gsimod.F90: Add commercial RO parameters + * guess_grids.F90: Add commercial RO parameters + * ncepnems_io.f90: Changes delz calculation to negative definite + * netcdfgfs_io.f90: Addresses GFSv16 bugzilla #1196 and delz sign bugfix + * read_files.f90: Addresses GFSv16 bugzilla #1196 + * read_fl_hdob.f90: Addresses GFSv16 bugzilla #1205 + * setupbend.f90: New QC and error inflation for commercial RO + * setuprad.f90: Bugfix for radiance diags + * stop1.f90: Addresses GFSv16 bugzilla #1196 + * write_incr.f90: Delz bugfix + + +USH CHANGES + +* ush/ + * build_all_cmake.sh: Addresses DA aspects of GFSv16 bugzillas + #216, #1218, and #1222 + * prune_4nco_global.sh: New script to address DA aspects of GFSv16 + bugzillas #216 and #1222 + + +UTIL CHANGES + +* util/Radiance_Monitor/nwprod/gdas_radmon.v3.0.0/scripts/exgdas_atmos_verfrad.sh: + chgrp Megha-Tropiques SAPHIR radiance diagnostic file to rsptrod since + this is a restricted data type, partially address DA aspects of GFSv16 + bugzilla #1221. + + +IMPLEMENTATION INSTRUCTIONS + +* The GFS DA v16.1 tag must be installed in conjunction with the entire + GFS v16.1 package. See release notes docs/Release_Notes.gfs.v16.1.0.txt + in the NOAA-EMC/global-workflow tag EMC-v16.1.0 tag for instructions. + + +PRE-IMPLEMENTATION TESTING REQUIREMENTS + +* Which production jobs should be tested as part of this implementation? + * The GFS DA package needs to be tested with the entire GFS suite. + +* Does this change require a 30-day evaluation? + * No. + + +DISSEMINATION INFORMATION + +* Where should this output be sent? + * No change from GFS v16.0 + +* Who are the users? + * No change from GFS v16.0 + +* Which output files should be transferred from PROD WCOSS to DEV WCOSS? + * No change from GFS v16.0 + +* Directory changes + * No change from GFS v16.0 + +* File changes + * No change from GFS v16.0 + + +HPSS ARCHIVE + +* No change from GFS v16.0 + + +JOB DEPENDENCIES & FLOW DIAGRAM + +* No change from GFSv16.0 + diff --git a/fix b/fix index daca7b586a..60d1b1fd3a 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit daca7b586aa66609cf5e1310c58695979617bc3f +Subproject commit 60d1b1fd3a73d9d8a5b7a743985911068eb5c6c2 diff --git a/scripts/exgdas_atmos_chgres_forenkf.sh b/scripts/exgdas_atmos_chgres_forenkf.sh index bb131c4276..416d58a2cc 100755 --- a/scripts/exgdas_atmos_chgres_forenkf.sh +++ b/scripts/exgdas_atmos_chgres_forenkf.sh @@ -48,7 +48,6 @@ export NMV=${NMV:-"/bin/mv"} export NLN=${NLN:-"/bin/ln -sf"} export CHGRP_CMD=${CHGRP_CMD:-"chgrp ${group_name:-rstprod}"} export NCLEN=${NCLEN:-$HOMEgfs/ush/getncdimlen} -export ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} # IAU DOIAU=${DOIAU:-"NO"} @@ -176,10 +175,7 @@ EOF . prep_step $APRUN_CHGRES $CHGRESNCEXEC chgres_nc_gauss0$FHR.nml - rc=$? - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 1 + export err=$?; err_chk fi done @@ -194,9 +190,7 @@ EOF . prep_step $APRUNCFP_CHGRES $DATA/mp_chgres.sh - export ERR=$? - export err=$ERR - $ERRSCRIPT || exit 3 + export err=$?; err_chk fi fi diff --git a/scripts/exgdas_enkf_ecen.sh b/scripts/exgdas_enkf_ecen.sh index 645baa3f7c..d15d568ebc 100755 --- a/scripts/exgdas_enkf_ecen.sh +++ b/scripts/exgdas_enkf_ecen.sh @@ -34,7 +34,6 @@ export CASE=${CASE:-384} ntiles=${ntiles:-6} # Utilities -ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} NCP=${NCP:-"/bin/cp -p"} NLN=${NLN:-"/bin/ln -sf"} NEMSIOGET=${NEMSIOGET:-${NWPROD}/exec/nemsio_get} @@ -164,11 +163,7 @@ if [ $DO_CALC_INCREMENT = "YES" ]; then $NCP $GETATMENSMEANEXEC $DATA $APRUN_ECEN ${DATA}/$(basename $GETATMENSMEANEXEC) $DATAPATH $ATMANLMEANNAME $ATMANLNAME $NMEM_ENKF - rc=$? - - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 2 + export err=$?; err_chk else # Link ensemble mean increment if [ $FHR -eq 6 ]; then @@ -188,11 +183,7 @@ else $NCP $GETATMENSMEANEXEC $DATA $APRUN_ECEN ${DATA}/$(basename $GETATMENSMEANEXEC) $DATAPATH $ATMINCMEANNAME $ATMINCNAME $NMEM_ENKF - rc=$? - - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 2 + export err=$?; err_chk # If available, link to ensemble mean guess. Otherwise, compute ensemble mean guess if [ -s $COMIN_GES_ENS/${GPREFIX}atmf00${FHR}.ensmean$GSUFFIX ]; then @@ -208,11 +199,7 @@ else $NCP $GETATMENSMEANEXEC $DATA $APRUN_ECEN ${DATA}/$(basename $GETATMENSMEANEXEC) $DATAPATH $ATMGESMEANNAME $ATMGESNAME $NMEM_ENKF - rc=$? - - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 2 + export err=$?; err_chk fi fi @@ -284,12 +271,7 @@ if [ $RECENTER_ENKF = "YES" ]; then EOF cat $chgresnml $APRUN_CHGRES ./chgres.x - rc=$? - - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 3 - + export err=$?; err_chk fi if [ $DO_CALC_INCREMENT = "YES" ]; then @@ -307,12 +289,7 @@ EOF $NCP $RECENATMEXEC $DATA $APRUN_ECEN ${DATA}/$(basename $RECENATMEXEC) $FILENAMEIN $FILENAME_MEANIN $FILENAME_MEANOUT $FILENAMEOUT $NMEM_ENKF - rc=$? - - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 2 - + export err=$?; err_chk else ################################################################################ # Recenter ensemble member atmospheric increments about hires analysis @@ -340,12 +317,7 @@ cat recenter.nml $NCP $RECENATMEXEC $DATA $APRUN_ECEN ${DATA}/$(basename $RECENATMEXEC) $FILENAMEIN $FILENAME_INCMEANIN $FILENAME_GSIDET $FILENAMEOUT $NMEM_ENKF $FILENAME_GESMEANIN - rc=$? - - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 2 - + export err=$?; err_chk fi fi @@ -388,11 +360,7 @@ EOF cat calc_increment.nml $APRUN_CALCINC ${DATA}/$(basename $CALCINCEXEC) - rc=$? - - export ERR=$rc - export err=$rc - $ERRSCRIPT || exit 4 + export err=$?; err_chk fi done # loop over analysis times in window diff --git a/scripts/exgdas_enkf_fcst.sh b/scripts/exgdas_enkf_fcst.sh index 45da2d4e78..2730055f19 100755 --- a/scripts/exgdas_enkf_fcst.sh +++ b/scripts/exgdas_enkf_fcst.sh @@ -34,8 +34,6 @@ export FIX_AM=${FIX_AM:-$FIX_DIR/fix_am} export NCP=${NCP:-"/bin/cp -p"} export NMV=${NMV:-"/bin/mv"} export NLN=${NLN:-"/bin/ln -sf"} -export ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} -export NDATE=${NDATE:-/$NWPROD/util/exec/ndate} # Scripts. FORECASTSH=${FORECASTSH:-$HOMEgfs/scripts/exglobal_forecast.sh} @@ -170,15 +168,13 @@ for imem in $(seq $ENSBEG $ENSEND); do export MEMBER=$imem export DATA=$DATATOP/$memchar if [ -d $DATA ]; then rm -rf $DATA; fi + mkdir -p $DATA $FORECASTSH ra=$? # Notify a member forecast failed and abort if [ $ra -ne 0 ]; then - msg="FATAL ERROR: forecast of member $cmem FAILED. Aborting job" - print $msg - export err=$ra - $ERRSCRIPT || exit 2 + err_exit "FATAL ERROR: forecast of member $cmem FAILED. Aborting job" fi ((rc+=ra)) @@ -226,9 +222,7 @@ cat $EFCSGRP ################################################################################ # If any members failed, error out -export ERR=$rc -export err=$ERR -$ERRSCRIPT || exit 2 +export err=$rc; err_chk ################################################################################ # Postprocessing diff --git a/scripts/exgdas_enkf_post.sh b/scripts/exgdas_enkf_post.sh index 32e709dc2f..2987e51d62 100755 --- a/scripts/exgdas_enkf_post.sh +++ b/scripts/exgdas_enkf_post.sh @@ -28,7 +28,6 @@ fi pwd=$(pwd) # Utilities -ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} NCP=${NCP:-"/bin/cp"} NLN=${NLN:-"/bin/ln -sf"} @@ -127,10 +126,7 @@ for fhr in $(seq $FHMIN $FHOUT $FHMAX); do ra=$? ((rc+=ra)) done - -export ERR=$rc -export err=$ERR -$ERRSCRIPT || exit 2 +export err=$rc; err_chk ################################################################################ # If smoothing on but no smoothing output, copy smoothed ensemble atmospheric files diff --git a/scripts/exgdas_enkf_select_obs.sh b/scripts/exgdas_enkf_select_obs.sh index 5298398e4f..a8a0cf1b22 100755 --- a/scripts/exgdas_enkf_select_obs.sh +++ b/scripts/exgdas_enkf_select_obs.sh @@ -29,7 +29,6 @@ pwd=$(pwd) # Utilities export NLN=${NLN:-"/bin/ln -sf"} -export ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} # Scripts. ANALYSISSH=${ANALYSISSH:-$HOMEgfs/scripts/exglobal_atmos_analysis.sh} @@ -115,10 +114,7 @@ export CHEM="$CHEM_INVOBS" # Execute GSI as a forward operator $ANALYSISSH - -export ERR=$? -export err=$ERR -$ERRSCRIPT || exit 2 +export err=$?; err_chk ################################################################################ # Postprocessing diff --git a/scripts/exgdas_enkf_sfc.sh b/scripts/exgdas_enkf_sfc.sh index f5ff74a824..4d927800d7 100755 --- a/scripts/exgdas_enkf_sfc.sh +++ b/scripts/exgdas_enkf_sfc.sh @@ -35,7 +35,6 @@ export CASE=${CASE:-384} ntiles=${ntiles:-6} # Utilities -ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} NCP=${NCP:-"/bin/cp -p"} NLN=${NLN:-"/bin/ln -sf"} NEMSIOGET=${NEMSIOGET:-${NWPROD}/exec/nemsio_get} @@ -163,10 +162,7 @@ if [ $DOIAU = "YES" ]; then done $CYCLESH - rc=$? - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 11 + export err=$?; err_chk done @@ -192,10 +188,7 @@ if [ $DOSFCANL_ENKF = "YES" ]; then done $CYCLESH - rc=$? - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 11 + export err=$?; err_chk done fi diff --git a/scripts/exgdas_enkf_update.sh b/scripts/exgdas_enkf_update.sh index 6e8186a55a..59b119fa13 100755 --- a/scripts/exgdas_enkf_update.sh +++ b/scripts/exgdas_enkf_update.sh @@ -30,7 +30,6 @@ pwd=$(pwd) # Utilities NCP=${NCP:-"/bin/cp -p"} NLN=${NLN:-"/bin/ln -sf"} -ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} NEMSIOGET=${NEMSIOGET:-$NWPROD/utils/exec/nemsio_get} NCLEN=${NCLEN:-$HOMEgfs/ush/getncdimlen} USE_CFP=${USE_CFP:-"NO"} @@ -245,9 +244,7 @@ if [ $USE_CFP = "YES" ]; then ncmd_max=$((ncmd < npe_node_max ? ncmd : npe_node_max)) APRUNCFP=$(eval echo $APRUNCFP) $APRUNCFP $DATA/mp_untar.sh - export ERR=$? - export err=$ERR - $ERRSCRIPT || exit 3 + export err=$?; err_chk fi fi @@ -378,11 +375,7 @@ export pgm=$ENKFEXEC $NCP $ENKFEXEC $DATA $APRUN_ENKF ${DATA}/$(basename $ENKFEXEC) 1>stdout 2>stderr -rc=$? - -export ERR=$rc -export err=$ERR -$ERRSCRIPT || exit 2 +export err=$?; err_chk # Cat runtime output files. cat stdout stderr > $COMOUT_ANL_ENS/$ENKFSTAT diff --git a/scripts/exglobal_atmos_analysis.sh b/scripts/exglobal_atmos_analysis.sh index ebca996f30..09656b8bd7 100755 --- a/scripts/exglobal_atmos_analysis.sh +++ b/scripts/exglobal_atmos_analysis.sh @@ -47,7 +47,6 @@ export NLN=${NLN:-"/bin/ln -sf"} export CHGRP_CMD=${CHGRP_CMD:-"chgrp ${group_name:-rstprod}"} export NEMSIOGET=${NEMSIOGET:-${NWPROD}/exec/nemsio_get} export NCLEN=${NCLEN:-$HOMEgfs/ush/getncdimlen} -export ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} COMPRESS=${COMPRESS:-gzip} UNCOMPRESS=${UNCOMPRESS:-gunzip} APRUNCFP=${APRUNCFP:-""} @@ -624,8 +623,7 @@ if [ $GENDIAG = "YES" ] ; then $NLN $DIAG_DIR/$pedir $pedir done else - echo "FATAL ERROR: lrun_subdirs must be true. lrun_subdirs=$lrun_subdirs" - $ERRSCRIPT || exit 2 + err_exit "FATAL ERROR: lrun_subdirs must be true. lrun_subdirs=$lrun_subdirs" fi fi @@ -716,9 +714,7 @@ EOFunzip ncmd_max=$((ncmd < npe_node_max ? ncmd : npe_node_max)) APRUNCFP_UNZIP=$(eval echo $APRUNCFP) $APRUNCFP_UNZIP $DATA/mp_unzip.sh - export ERR=$? - export err=$ERR - $ERRSCRIPT || exit 3 + export err=$?; err_chk fi fi fi # if [ $USE_RADSTAT = "YES" ] @@ -753,7 +749,7 @@ cat > gsiparm.anl << EOF iguess=-1, tzr_qc=$TZR_QC, oneobtest=.false.,retrieval=.false.,l_foto=.false., - use_pbl=.false.,use_compress=.true.,nsig_ext=12,gpstop=50., + use_pbl=.false.,use_compress=.true.,nsig_ext=12,gpstop=50.,commgpstop=45.,commgpserrinf=1.0, use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},sfcnst_comb=.true., use_readin_anl_sfcmask=${USE_READIN_ANL_SFCMASK}, lrun_subdirs=$lrun_subdirs, @@ -958,11 +954,7 @@ export pgm=$GSIEXEC $NCP $GSIEXEC $DATA $APRUN_GSI ${DATA}/$(basename $GSIEXEC) 1>&1 2>&2 -rc=$? - -export ERR=$rc -export err=$ERR -$ERRSCRIPT || exit 4 +export err=$?; err_chk ############################################################## @@ -970,11 +962,7 @@ $ERRSCRIPT || exit 4 # here before releasing FV3 forecast if [ $DO_CALC_INCREMENT = "YES" ]; then $CALCINCPY - rc=$? - - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 5 + export err=$?; err_chk fi ############################################################## @@ -1030,10 +1018,7 @@ if [ $DOGCYCLE = "YES" ]; then export MAX_TASKS_CY=$ntiles $CYCLESH - rc=$? - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 11 + export err=$?; err_chk fi # update surface restarts at middle of window for n in $(seq 1 $ntiles); do @@ -1048,11 +1033,7 @@ if [ $DOGCYCLE = "YES" ]; then export MAX_TASKS_CY=$ntiles $CYCLESH - rc=$? - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 11 - + export err=$?; err_chk fi @@ -1072,6 +1053,7 @@ if [ $RUN_SELECT = "YES" ]; then echo $(date) START tar obs_input >&2 [[ -s obsinput.tar ]] && rm obsinput.tar $NLN $SELECT_OBS obsinput.tar + ${CHGRP_CMD} obs_input.* tar -cvf obsinput.tar obs_input.* chmod 750 $SELECT_OBS ${CHGRP_CMD} $SELECT_OBS diff --git a/scripts/exglobal_atmos_analysis_calc.sh b/scripts/exglobal_atmos_analysis_calc.sh index 7204af81f0..13e53874fa 100755 --- a/scripts/exglobal_atmos_analysis_calc.sh +++ b/scripts/exglobal_atmos_analysis_calc.sh @@ -48,7 +48,6 @@ export NLN=${NLN:-"/bin/ln -sf"} export CHGRP_CMD=${CHGRP_CMD:-"chgrp ${group_name:-rstprod}"} export NEMSIOGET=${NEMSIOGET:-${NWPROD}/exec/nemsio_get} export NCLEN=${NCLEN:-$HOMEgfs/ush/getncdimlen} -export ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} COMPRESS=${COMPRESS:-gzip} UNCOMPRESS=${UNCOMPRESS:-gunzip} APRUNCFP=${APRUNCFP:-""} @@ -186,11 +185,7 @@ if [ $DO_CALC_ANALYSIS == "YES" ]; then fi $CALCANLPY - rc=$? - - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 3 + export err=$?; err_chk else echo "Neither increment nor analysis are generated by external utils" fi @@ -202,10 +197,7 @@ if [ $DOGAUSFCANL = "YES" ]; then export OMP_NUM_THREADS_SFC=$NTHREADS_GAUSFCANL $GAUSFCANLSH - rc=$? - export ERR=$rc - export err=$ERR - $ERRSCRIPT || exit 12 + export err=$?; err_chk fi echo "$CDUMP $CDATE atmanl and sfcanl done at `date`" > $COMOUT/${APREFIX}loganl.txt diff --git a/scripts/exglobal_diag.sh b/scripts/exglobal_diag.sh index 1afc49f246..eefc895b86 100755 --- a/scripts/exglobal_diag.sh +++ b/scripts/exglobal_diag.sh @@ -48,7 +48,6 @@ export CHGRP_CMD=${CHGRP_CMD:-"chgrp ${group_name:-rstprod}"} export NEMSIOGET=${NEMSIOGET:-${NWPROD}/exec/nemsio_get} export NCLEN=${NCLEN:-$HOMEgfs/ush/getncdimlen} export CATEXEC=${CATEXEC:-$HOMEgfs/exec/ncdiag_cat.x} -export ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} COMPRESS=${COMPRESS:-gzip} UNCOMPRESS=${UNCOMPRESS:-gunzip} APRUNCFP=${APRUNCFP:-""} @@ -108,10 +107,7 @@ if [ $GENDIAG = "YES" ] ; then $NLN $pe $DATA/$pedir done else - echo "lrun_subdirs must be true; exit with error" - export ERR=$? - export err=$ERR - $ERRSCRIPT || exit 2 + err_exit "***FATAL ERROR*** lrun_subdirs must be true. Abort job" fi # Set up lists and variables for various types of diagnostic files. @@ -245,12 +241,16 @@ EOFdiag ncmd_max=$((ncmd < npe_node_max ? ncmd : npe_node_max)) APRUNCFP_DIAG=$(eval echo $APRUNCFP) $APRUNCFP_DIAG $DATA/mp_diag.sh - export ERR=$? - export err=$ERR - $ERRSCRIPT || exit 3 + export err=$?; err_chk fi fi + # Restrict diagnostic files containing rstprod data + rlist="conv_gps conv_ps conv_pw conv_q conv_sst conv_t conv_uv saphir" + for rtype in $rlist; do + ${CHGRP_CMD} *${rtype}* + done + # If requested, create diagnostic file tarballs if [ $DIAG_TARBALL = "YES" ]; then echo $(date) START tar diagnostic files >&2 @@ -262,9 +262,7 @@ EOFdiag fi if [ ${numfile[n]} -gt 0 ]; then tar $TAROPTS ${diagfile[n]} $(cat ${diaglist[n]}) - export ERR=$? - export err=$ERR - $ERRSCRIPT || exit 4 + export err=$?; err_chk fi done diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index f245afdf3f..58ceff78a6 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -207,6 +207,7 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & endif ! pressure at interfaces do k=1,nlevs+1 + ! k=1 in ak,bk is model top pressi(:,k) = 0.01_r_kind*ak(nlevs-k+2)+bk(nlevs-k+2)*psg if (nanal .eq. 1 .and. iope==0) print *,'netcdf, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) enddo @@ -624,6 +625,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, endif ! pressure at interfaces do k=1,nlevs+1 + ! k=1 in ak,bk is model top pressi(:,k) = 0.01_r_kind*ak(nlevs-k+2)+bk(nlevs-k+2)*psg if (nanal .eq. 1) print *,'netcdf, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) enddo @@ -1146,6 +1148,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_ call stop2(29) endif do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom ak(nlevs-k+2) = 0.01_r_kind*values_1d(k) enddo call read_attribute(dsfg, 'bk', values_1d,errcode=iret) @@ -1154,6 +1157,7 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_ call stop2(29) endif do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom bk(nlevs-k+2) = values_1d(k) enddo @@ -1509,16 +1513,18 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_ endif vg = values_1d + (vg*100_r_kind) ! analysis ps (values_1d is background ps) do k=lev_pe1(iope), lev_pe2(iope) - krev = nlevs-k+1 + krev = nlevs-k+1 ! k=1 is model top ki = k - lev_pe1(iope) + 1 ug=(rd/grav)*reshape(tv_anal(:,:,ki),(/nlons*nlats/)) ! ps in Pa here, need to multiply ak by 100. - ug=ug*log((100_r_kind*ak(krev)+bk(krev)*vg)/(100_r_kind*ak(krev+1)+bk(krev+1)*vg)) + ! ug is hydrostatic analysis delz (should be negative) + ! (note that ak,bk are already reversed to go from bottom to top) + ug=ug*log((100_r_kind*ak(krev+1)+bk(krev+1)*vg)/(100_r_kind*ak(krev)+bk(krev)*vg)) ! ug is hydrostatic analysis delz inferred from analysis ps,Tv - ! delzb is hydrostatic background delz inferred from background ps,Tv + ! delzb is (negative) hydrostatic background delz inferred from background ps,Tv delzb=(rd/grav)*reshape(tv_bg(:,:,ki),(/nlons*nlats/)) - delzb=delzb*log((100_r_kind*ak(krev)+bk(krev)*values_1d)/(100_r_kind*ak(krev+1)+bk(krev+1)*values_1d)) - ug3d(:,:,ki)=values_3d(:,:,ki) + reshape(delzb-ug,(/nlons,nlats/)) + delzb=delzb*log((100_r_kind*ak(krev+1)+bk(krev+1)*values_1d)/(100_r_kind*ak(krev)+bk(krev)*values_1d)) + ug3d(:,:,ki)=values_3d(:,:,ki) + reshape(ug-delzb,(/nlons,nlats/)) end do if (has_attr(dsfg, 'nbits', 'delz') .and. .not. nocompress) then call read_attribute(dsfg, 'nbits', nbits, 'delz') @@ -1881,7 +1887,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n integer :: ps_ind, pst_ind, nbits integer :: ql_ind, qi_ind, qr_ind, qs_ind, qg_ind - integer k,nt,ierr,iunitsig,nb,i,ne,nanal + integer k,krev,nt,ierr,iunitsig,nb,i,ne,nanal logical :: nocompress @@ -1969,6 +1975,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n call stop2(29) endif do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom ak(nlevs-k+2) = 0.01_r_kind*values_1d(k) enddo call read_attribute(dsfg, 'bk', values_1d,errcode=iret) @@ -1977,6 +1984,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n call stop2(29) endif do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom bk(nlevs-k+2) = values_1d(k) enddo else @@ -2475,8 +2483,9 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n if (hasfield) then call nemsio_readrecv(gfilein,'pres','sfc',1,nems_wrk2,iret=iret) delzb=(rd/grav)*nems_wrk - ! ps in Pa here, need to multiply ak by 100. - delzb=delzb*log((100_r_kind*ak(k)+bk(k)*nems_wrk2)/(100_r_kind*ak(k+1)+bk(k+1)*nems_wrk2)) + ! ps in Pa here, need to multiply ak by 100. k is bottom to top, + ! calculate delzb so it is negative. + delzb=delzb*log((100_r_kind*ak(k+1)+bk(k+1)*nems_wrk2)/(100_r_kind*ak(k)+bk(k)*nems_wrk2)) endif ! convert Tv back to T nems_wrk = ug/(1. + fv*vg) @@ -2508,14 +2517,16 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n vg = nems_wrk2 + vg ug=(rd/grav)*ug ! ug is analysis Tv ! ps in Pa here, need to multiply ak by 100. - ug=ug*log((100_r_kind*ak(k)+bk(k)*vg)/(100_r_kind*ak(k+1)+bk(k+1)*vg)) - ug=ug-delzb + ! k is bottom to top, calculate delz so it is negative + ug=ug*log((100_r_kind*ak(k+1)+bk(k+1)*vg)/(100_r_kind*ak(k)+bk(k)*vg)) + ug=ug-delzb ! analysis - background call nemsio_readrecv(gfilein,'delz','mid layer',k,nems_wrk,iret=iret) if (iret/=0) then write(6,*)'gridio/writegriddata: gfs model: problem with nemsio_readrecv(delz), iret=',iret call stop2(23) endif - if (sum(nems_wrk) < 0.0_r_kind) ug = ug * -1.0_r_kind + ! flip sign of delz increment if background is positive. + if (sum(nems_wrk) > 0.0_r_kind) ug = ug * -1.0_r_kind nems_wrk = nems_wrk + ug call nemsio_writerecv(gfileout,'delz','mid layer',k,nems_wrk,iret=iret) if (iret/=0) then @@ -2876,17 +2887,20 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n endif vg = values_1d + vg*100_r_kind ! analysis ps (values_1d is background ps) do k=1,nlevs - ug=(rd/grav)*reshape(tv_anal(:,:,nlevs-k+1),(/nlons*nlats/)) + krev = nlevs-k+1 ! k=1 is model top + ug=(rd/grav)*reshape(tv_anal(:,:,k),(/nlons*nlats/)) ! ps in Pa here, need to multiply ak by 100. - ug=ug*log((100_r_kind*ak(k)+bk(k)*vg)/(100_r_kind*ak(k+1)+bk(k+1)*vg)) + ! ug is analysis delz, calculate so it is negative + ! (note that ak,bk are already reversed to go from bottom to top) + ug=ug*log((100_r_kind*ak(krev+1)+bk(krev+1)*vg)/(100_r_kind*ak(krev)+bk(krev)*vg)) ! ug is hydrostatic analysis delz inferred from analysis ps,Tv ! delzb is hydrostatic background delz inferred from background ps,Tv - delzb=(rd/grav)*reshape(tv_bg(:,:,nlevs-k+1),(/nlons*nlats/)) - delzb=delzb*log((100_r_kind*ak(k)+bk(k)*values_1d)/(100_r_kind*ak(k+1)+bk(k+1)*values_1d)) - ug3d(:,:,nlevs-k+1)=values_3d(:,:,nlevs-k+1) +& - reshape(delzb-ug,(/nlons,nlats/)) + ! calculate so it is negative + delzb=(rd/grav)*reshape(tv_bg(:,:,k),(/nlons*nlats/)) + delzb=delzb*log((100_r_kind*ak(krev+1)+bk(krev+1)*values_1d)/(100_r_kind*ak(krev)+bk(krev)*values_1d)) + ug3d(:,:,k)=values_3d(:,:,k) +& + reshape(ug-delzb,(/nlons,nlats/)) enddo - !print *,'min/max delz',minval(values_3d),maxval(values_3d),& ! minval(ug3d),maxval(ug3d) if (has_attr(dsfg, 'nbits', 'delz') .and. .not. nocompress) then call read_attribute(dsfg, 'nbits', nbits, 'delz') @@ -3409,6 +3423,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, call stop2(29) endif do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom ak(nlevs-k+2) = 0.01_r_kind*values_1d(k) enddo call read_attribute(dsfg, 'bk', values_1d,errcode=iret) @@ -3417,6 +3432,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, call stop2(29) endif do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom bk(nlevs-k+2) = values_1d(k) enddo @@ -3571,15 +3587,17 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, psges = reshape(values_2d,(/nlons*nlats/)) vg = psges + (psinc*100_r_kind) do k=1,nlevs - krev = nlevs-k+1 + krev = nlevs-k+1 ! k=1 is model top ug=(rd/grav)*reshape(tvanl(:,:,k),(/nlons*nlats/)) ! ps in Pa here, need to multiply ak by 100. - ug=ug*log((100_r_kind*ak(krev)+bk(krev)*vg)/(100_r_kind*ak(krev+1)+bk(krev+1)*vg)) + ! calculate ug (analysis delz) so it is negative. + ! (note that ak,bk are already reversed to go from bottom to top) + ug=ug*log((100_r_kind*ak(krev+1)+bk(krev+1)*vg)/(100_r_kind*ak(krev)+bk(krev)*vg)) ! ug is hydrostatic analysis delz inferred from analysis ps,Tv ! delzb is hydrostatic background delz inferred from background ps,Tv delzb=(rd/grav)*reshape(tv(:,:,k),(/nlons*nlats/)) - delzb=delzb*log((100_r_kind*ak(krev)+bk(krev)*psges)/(100_r_kind*ak(krev+1)+bk(krev+1)*psges)) - inc3d(:,:,k)=reshape(delzb-ug,(/nlons,nlats/)) + delzb=delzb*log((100_r_kind*ak(krev+1)+bk(krev+1)*psges)/(100_r_kind*ak(krev)+bk(krev)*psges)) + inc3d(:,:,k)=reshape(ug-delzb,(/nlons,nlats/)) end do end if do j=1,nlats @@ -3841,6 +3859,7 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate call stop2(29) endif do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom ak(nlevs-k+2) = 0.01_r_kind*values_1d(k) enddo call read_attribute(dsfg, 'bk', values_1d,errcode=iret) @@ -3849,6 +3868,7 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate call stop2(29) endif do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom bk(nlevs-k+2) = values_1d(k) enddo @@ -4016,12 +4036,15 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate ki = k - lev_pe1(iope) + 1 ug=(rd/grav)*reshape(tvanl(:,:,ki),(/nlons*nlats/)) ! ps in Pa here, need to multiply ak by 100. - ug=ug*log((100_r_kind*ak(krev)+bk(krev)*vg)/(100_r_kind*ak(krev+1)+bk(krev+1)*vg)) + ! calculate analysis delz so it is negative. + ! (note that ak,bk are already reversed to go from bottom to top) + ug=ug*log((100_r_kind*ak(krev+1)+bk(krev+1)*vg)/(100_r_kind*ak(krev)+bk(krev)*vg)) ! ug is hydrostatic analysis delz inferred from analysis ps,Tv ! delzb is hydrostatic background delz inferred from background ps,Tv + ! calculate delzb so it is negative delzb=(rd/grav)*reshape(tv(:,:,ki),(/nlons*nlats/)) - delzb=delzb*log((100_r_kind*ak(krev)+bk(krev)*psges)/(100_r_kind*ak(krev+1)+bk(krev+1)*psges)) - inc3d(:,:,ki)=reshape(delzb-ug,(/nlons,nlats/)) + delzb=delzb*log((100_r_kind*ak(krev+1)+bk(krev+1)*psges)/(100_r_kind*ak(krev)+bk(krev)*psges)) + inc3d(:,:,ki)=reshape(ug-delzb,(/nlons,nlats/)) end do end if do j=1,nlats diff --git a/src/fv3gfs_ncio/module_fv3gfs_ncio.f90 b/src/fv3gfs_ncio/module_fv3gfs_ncio.f90 index e7a38b63b0..51b901bf88 100644 --- a/src/fv3gfs_ncio/module_fv3gfs_ncio.f90 +++ b/src/fv3gfs_ncio/module_fv3gfs_ncio.f90 @@ -111,7 +111,7 @@ subroutine nccheck(status,halt,fname) implicit none integer, intent (in) :: status logical, intent(in), optional :: halt - character(len=500), intent(in), optional :: fname + character(len=*), intent(in), optional :: fname logical stopit if (present(halt)) then stopit = halt diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90 index 8baf06067a..ac00db4c41 100644 --- a/src/gsi/cplr_gfs_ensmod.f90 +++ b/src/gsi/cplr_gfs_ensmod.f90 @@ -914,7 +914,7 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig character(len=*), intent(in ) :: filename ! Declare local variables - integer(i_kind) i,ii,j,jj,k,lonb,latb,levs,kr + integer(i_kind) i,ii,j,jj,k,lonb,latb,levs,kr,ierror integer(i_kind) k2,k3,k3u,k3v,k3t,k3q,k3cw,k3oz,kf character(len=120) :: myname_ = 'parallel_read_gfsnc_state_' real(r_single),allocatable,dimension(:,:,:) :: rwork3d1, rwork3d2 @@ -927,7 +927,11 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig type(Dimension) :: ncdim - atmges = open_dataset(filename) + atmges = open_dataset(filename,errcode=ierror) + if (ierror /=0) then + write(6,*)' PARALLEL_READ_GFSNC_STATE: ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS' + call stop2(999) + endif ! get dimension sizes ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len diff --git a/src/gsi/gesinfo.F90 b/src/gsi/gesinfo.F90 index 48757c2135..5d86c04dc6 100644 --- a/src/gsi/gesinfo.F90 +++ b/src/gsi/gesinfo.F90 @@ -337,8 +337,16 @@ subroutine gesinfo else ! use_gfs_ncio and get this information write(sfilename,'("sfcf",i2.2)')nhr_assimilation ! open the netCDF file - atmges = open_dataset(filename) - sfcges = open_dataset(sfilename) + atmges = open_dataset(filename,errcode=iret) + if (iret /=0) then + write(6,*)'GESINFO: ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS' + call stop2(99) + endif + sfcges = open_dataset(sfilename,errcode=iret) + if (iret /=0) then + write(6,*)'GESINFO: ***ERROR*** ',trim(sfilename),' NOT AVAILABLE: PROGRAM STOPS' + call stop2(99) + endif ! get dimension sizes ncdim = get_dim(atmges, 'grid_xt'); gfshead%lonb = ncdim%len ncdim = get_dim(atmges, 'grid_yt'); gfshead%latb = ncdim%len diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index e1cc923f18..78bebc5140 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -124,7 +124,7 @@ module gsimod jcap_gfs,nlat_gfs,nlon_gfs,jcap_cut,wrf_mass_hybridcord,use_gfs_ncio,write_fv3_incr,& use_fv3_aero,grid_type_fv3_regional use gridmod,only: l_reg_update_hydro_delz - use guess_grids, only: ifact10,sfcmod_gfs,sfcmod_mm5,use_compress,nsig_ext,gpstop + use guess_grids, only: ifact10,sfcmod_gfs,sfcmod_mm5,use_compress,nsig_ext,gpstop,commgpstop,commgpserrinf use gsi_io, only: init_io,lendian_in,verbose,print_obs_para use regional_io_mod, only: regional_io_class use wrf_params_mod, only: update_pint, preserve_restart_date @@ -572,6 +572,8 @@ module gsimod ! use_compress - option to turn on the use of compressibility factors in geopotential heights ! nsig_ext - number of layers above the model top which are necessary to compute the bending angle for gpsro ! gpstop - maximum height for gpsro data assimilation. Reject anything above this height. +! commgpstop -Reject commercial ro above this height. Logic in setupbend assumes commgpstop <= gpstop. +! commgpserrinf - optional error inflation factor for commercial gpsro data ! use_gfs_nemsio - option to use nemsio to read global model NEMS/GFS first guess ! use_gfs_ncio - option to use netCDF to read global model FV3-GFS first guess ! use_fv3_aero - option to use FV3-Chem vs NGAC for global aerosol analysis @@ -664,7 +666,7 @@ module gsimod diag_rad,diag_pcp,diag_conv,diag_ozone,diag_aero,diag_co,diag_light,diag_radardbz,iguess, & write_diag,reduce_diag, & oneobtest,sfcmodel,dtbduv_on,ifact10,l_foto,offtime_data,& - use_pbl,use_compress,nsig_ext,gpstop,& + use_pbl,use_compress,nsig_ext,gpstop,commgpstop, commgpserrinf, & perturb_obs,perturb_fact,oberror_tune,preserve_restart_date, & crtm_coeffs_path,berror_stats,tcp_posmatch,tcp_box, & newpc4pred,adp_anglebc,angord,passive_bc,use_edges,emiss_bc,upd_pred,reset_bad_radbc,& diff --git a/src/gsi/guess_grids.F90 b/src/gsi/guess_grids.F90 index 2d3d55637f..27d9c7a466 100644 --- a/src/gsi/guess_grids.F90 +++ b/src/gsi/guess_grids.F90 @@ -144,7 +144,7 @@ module guess_grids public :: geom_hgti,geom_hgti_bg public :: wgt_lcbas public :: ges_qsat - public :: use_compress,nsig_ext,gpstop + public :: use_compress,nsig_ext,gpstop,commgpstop,commgpserrinf public :: ges_tsen1,ges_q1 public :: ntguesaer,ifileaer,nfldaer,hrdifaer ! variables for external aerosol files @@ -230,7 +230,8 @@ module guess_grids real(r_kind):: gpstop=30.0_r_kind ! maximum gpsro height used in km ! geometric height for ref, impact height for bnd - + real(r_kind):: commgpstop=30.0_r_kind + real(r_kind):: commgpserrinf=1.0_r_kind ! error inflation factor for commercial gnssro real(r_kind):: ges_psfcavg ! average guess surface pressure real(r_kind),allocatable,dimension(:):: ges_prslavg ! average guess pressure profile diff --git a/src/gsi/ncepnems_io.f90 b/src/gsi/ncepnems_io.f90 index 0525b83bf9..b4c29087d7 100755 --- a/src/gsi/ncepnems_io.f90 +++ b/src/gsi/ncepnems_io.f90 @@ -2310,12 +2310,13 @@ subroutine write_fv3atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) ! Calculate delz increment for UPP if (lupp) then do k=1,grd%nsig - sub_dzb(:,:,k) = ges_geopi(:,:,k+1,ibin) - ges_geopi(:,:,k,ibin) + ! k goes from bottom to top, so k - k+1 is negative delz + sub_dzb(:,:,k) = ges_geopi(:,:,k,ibin) - ges_geopi(:,:,k+1,ibin) enddo if ((.not. lwrite4danl) .or. ibin == 1) call load_geop_hgt do k=1,grd%nsig - sub_dza(:,:,k) = geop_hgti(:,:,k+1,ibin) - geop_hgti(:,:,k,ibin) + sub_dza(:,:,k) = geop_hgti(:,:,k,ibin) - geop_hgti(:,:,k+1,ibin) enddo sub_dza = sub_dza - sub_dzb !sub_dza is increment @@ -2800,10 +2801,10 @@ subroutine write_fv3atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then - work1 = -one * work1 ! Flip sign, FV3 is top to bottom call nemsio_readrecv(gfile,'delz','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','read',istop,iret) - if (sum(rwork1d) < zero) work1 = -one * work1 !Flip sign, FV3 is top to bottom + ! if delz in background is positive, flip sign of increment + if (sum(rwork1d) > zero) work1 = -one * work1 if(diff_res)then grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) do kk=1,grd%iglobal @@ -3149,12 +3150,13 @@ subroutine write_atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin,gfschem_bundle ! Calculate delz increment for UPP if (lupp) then do k=1,grd%nsig - sub_dzb(:,:,k) = ges_geopi(:,:,k+1,ibin) - ges_geopi(:,:,k,ibin) + ! k goes from bottom to top, so k - k+1 is negative delz + sub_dzb(:,:,k) = ges_geopi(:,:,k,ibin) - ges_geopi(:,:,k+1,ibin) enddo if ((.not. lwrite4danl) .or. ibin == 1) call load_geop_hgt do k=1,grd%nsig - sub_dza(:,:,k) = geop_hgti(:,:,k+1,ibin) - geop_hgti(:,:,k,ibin) + sub_dza(:,:,k) = geop_hgti(:,:,k,ibin) - geop_hgti(:,:,k+1,ibin) enddo sub_dza = sub_dza - sub_dzb !sub_dza is increment @@ -3531,10 +3533,10 @@ subroutine write_atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin,gfschem_bundle work1,grd%ijn,grd%displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype == mype_out) then - work1 = -one * work1 ! Flip sign, FV3 is top to bottom call nemsio_readrecv(gfile,'delz','mid layer',k,rwork1d,iret=iret) if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','read',istop,iret) - if (sum(rwork1d) < zero) work1 = -one * work1 ! Flip sign, FV3 is top to bottom + ! if delz in background is positive, flip sign of increment + if (sum(rwork1d) > zero) work1 = -one * work1 ! Flip sign, FV3 is top to bottom if(diff_res)then grid_b=reshape(rwork1d,(/size(grid_b,1),size(grid_b,2)/)) do kk=1,grd%iglobal diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index eb90b7d142..dd3f991da5 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -882,7 +882,7 @@ subroutine read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, & real(r_single), dimension(nlat_sfc,nlon_sfc,nfldsfc) :: xt character(len=24) :: filename character(len=120) :: my_name = 'READ_SFCNST' - integer(i_kind) :: i,j,it,n,nsfc + integer(i_kind) :: i,j,it,n,nsfc,iret integer(i_kind) :: lonb, latb real(r_single),allocatable, dimension(:) :: fhour real(r_single), allocatable, dimension(:,:) :: work,outtmp @@ -896,7 +896,11 @@ subroutine read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, & 200 format('sfcf',i2.2) ! open the netCDF file - sfcges = open_dataset(filename) + sfcges = open_dataset(filename,errcode=iret) + if (iret/=0) then + write(6,*) trim(my_name),': ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS, later' + endif + ! get dimension sizes ncdim = get_dim(sfcges, 'grid_xt'); lonb = ncdim%len ncdim = get_dim(sfcges, 'grid_yt'); latb = ncdim%len @@ -1212,7 +1216,7 @@ subroutine read_sfc_anl_(isli_anl) ! Declare local variables character(len=24) :: filename character(len=120) :: my_name = 'READ_GFSNCSFC_ANL' - integer(i_kind) :: i,j + integer(i_kind) :: i,j,iret integer(i_kind) :: lonb, latb real(r_single),allocatable, dimension(:) :: fhour real(r_single), allocatable, dimension(:,:) :: work,outtmp @@ -1223,7 +1227,12 @@ subroutine read_sfc_anl_(isli_anl) filename='sfcf06_anlgrid' ! open the netCDF file - sfcges = open_dataset(filename) + sfcges = open_dataset(filename,errcode=iret) + if (iret/=0) then + write(6,*) trim(my_name),': ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS' + call stop2(999) + endif + ! get dimension sizes ncdim = get_dim(sfcges, 'grid_xt'); lonb = ncdim%len ncdim = get_dim(sfcges, 'grid_yt'); latb = ncdim%len @@ -1774,14 +1783,14 @@ subroutine write_atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) endif ! if ( mype == mype_out ) - ! Calculate delz increment for UPP + ! compute delz (so that delz < 0 as model expects) do k=1,grd%nsig - sub_dzb(:,:,k) = ges_geopi(:,:,k+1,ibin) - ges_geopi(:,:,k,ibin) + sub_dzb(:,:,k) = ges_geopi(:,:,k,ibin) - ges_geopi(:,:,k+1,ibin) enddo if ((.not. lwrite4danl) .or. ibin == 1) call load_geop_hgt do k=1,grd%nsig - sub_dza(:,:,k) = geop_hgti(:,:,k+1,ibin) - geop_hgti(:,:,k,ibin) + sub_dza(:,:,k) = geop_hgti(:,:,k,ibin) - geop_hgti(:,:,k+1,ibin) enddo sub_dza = sub_dza - sub_dzb !sub_dza is increment @@ -2201,7 +2210,7 @@ subroutine write_atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) values_3d(:,:,kr) = grid_b else call load_grid(work1,grid) - values_3d(:,:,kr) = values_3d(:,:,kr) - grid + values_3d(:,:,kr) = values_3d(:,:,kr) + grid end if end if endif @@ -2378,6 +2387,7 @@ subroutine write_sfc_ (filename,mype_sfc,dsfct) ! Read surface guess file ! open the netCDF file sfcges = open_dataset(fname_ges) + ! get dimension sizes ncdim = get_dim(sfcges, 'grid_xt'); lonb = ncdim%len ncdim = get_dim(sfcges, 'grid_yt'); latb = ncdim%len diff --git a/src/gsi/read_files.f90 b/src/gsi/read_files.f90 index 5fa4119c7b..6f1e7a1488 100644 --- a/src/gsi/read_files.f90 +++ b/src/gsi/read_files.f90 @@ -82,7 +82,7 @@ subroutine read_files(mype) ifilesig,ifilesfc,ifilenst,hrdifsig,hrdifsfc,hrdifnst,create_gesfinfo use guess_grids, only: hrdifsig_all,hrdifsfc_all,hrdifnst_all use guess_grids, only: nfldaer, ntguesaer, ifileaer, hrdifaer, hrdifaer_all !for aerosol - use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen,nhr_assimilation + use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen,nhr_assimilation,nhr_obsbin use hybrid_ensemble_parameters, only: ntlevs_ens use gridmod, only: nlat_sfc,nlon_sfc,lpl_gfs,dx_gfs,use_gfs_nemsio,sfcnst_comb,use_gfs_ncio use constants, only: zero,r60inv,r60,r3600,i_missing @@ -118,9 +118,10 @@ subroutine read_files(mype) ! Declare local variables logical(4) fexist + logical:: present character(6) filename integer(i_kind) i,j,iwan,npem1,iret - integer(i_kind) nhr_half + integer(i_kind) nhr_half,ihr integer(i_kind) iamana(4) ! changed to 4 from 3 for aer files integer(i_kind) nminanl,nmings,nming2,ndiff integer(i_kind),dimension(4):: idateg @@ -258,7 +259,9 @@ subroutine read_files(mype) idateg=sigatm_head%idate call sigio_sclose(lunatm,iret) else if (use_gfs_ncio) then - atmges = open_dataset(filename) + atmges = open_dataset(filename,errcode=iret) + if (iret /=0 .and. mype==0) & + write(6,*)'READ_FILES: ***WARNING*** problem reading atm file ',trim(filename),iret idate6 = get_idate_from_time_units(atmges) call read_vardata(atmges, 'time', fhour) hourg4 = float(nint(fhour(1))) ! going to make this nearest integer for now @@ -332,7 +335,9 @@ subroutine read_files(mype) call sfcio_sclose(lunsfc,iret) if(i == 1 .and. print_verbose)write(6,*)' READ_FILES: in sfcio sfc_head%lpl = ', sfc_head%lpl else if (use_gfs_ncio) then - sfcges = open_dataset(filename) + sfcges = open_dataset(filename,errcode=iret) + if (iret /=0 .and. mype==0) & + write(6,*)'READ_FILES: ***WARNING*** problem reading sfc file ',trim(filename),iret ncdim = get_dim(sfcges, 'grid_xt'); sfc_head%lonb = ncdim%len ncdim = get_dim(sfcges, 'grid_yt'); sfc_head%latb = ncdim%len idate6 = get_idate_from_time_units(sfcges) @@ -614,7 +619,21 @@ subroutine read_files(mype) call stop2(99) endif if (l4densvar .and. nfldsig/=ntlevs_ens) then - write(6,*)'READ_FILES: ***ERROR*** insufficient atm fcst for 4densvar: PROGRAM STOPS' + if (mype==0) then + write(6,*)'READ_FILES: ***ERROR*** insufficient atm fcst for 4densvar: PROGRAM STOPS' + do i=1,ntlevs_ens + ihr=nhr_obsbin*(i-1)+nhr_half + present=.false. + do j=1,nfldsig + if (ihr == ifilesig(j)) present=.true. + end do + if (.not.present) then + write(filename,'(''sigf'',i2.2)')ihr + write(6,*)'READ_FILES: ***ERROR*** file ',trim(filename),' missing: PROGRAM STOPS' + endif + end do + endif + call mpi_barrier(mpi_comm_world,ierror) call stop2(99) endif @@ -632,7 +651,21 @@ subroutine read_files(mype) call stop2(99) endif if (l4densvar .and. nfldsfc/=ntlevs_ens) then - write(6,*)'READ_FILES: ***ERROR*** insufficient sfc fcst for 4densvar: PROGRAM STOPS' + if (mype==0) then + write(6,*)'READ_FILES: ***ERROR*** insufficient sfc fcst for 4densvar: PROGRAM STOPS' + do i=1,ntlevs_ens + ihr=nhr_obsbin*(i-1)+nhr_half + present=.false. + do j=1,nfldsfc + if (ihr == ifilesfc(j)) present=.true. + end do + if (.not.present) then + write(filename,'(''sfcf'',i2.2)')ihr + write(6,*)'READ_FILES: ***ERROR*** file ',trim(filename),' missing: PROGRAM STOPS' + endif + end do + endif + call mpi_barrier(mpi_comm_world,ierror) call stop2(99) endif diff --git a/src/gsi/read_fl_hdob.f90 b/src/gsi/read_fl_hdob.f90 index e600cb7bae..c7dc95f612 100644 --- a/src/gsi/read_fl_hdob.f90 +++ b/src/gsi/read_fl_hdob.f90 @@ -769,8 +769,6 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si tdob_calc = tob*(one-tob*log(rhob/100)) ! for comparison qob = rhob*qsat endif -! write(4000,1004) nread,pob_mb,tob,tdob,qob,qsat,rhob,q_qm,usage -!1004 format(i6,6(1x,e20.12),1x,i5,1x,f5.0) ! Get observation error from error table if (njqc) then ppb = max(zero,min(pob_mb,r2000)) @@ -1043,8 +1041,6 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si cdata_all(26,iout)=ran01dom()*perturb_fact ! u perturbation cdata_all(27,iout)=ran01dom()*perturb_fact ! v perturbation endif - write(3000,1003) nread,pob_mb,uob,vob,qcm,usage -1003 format(i12,3e25.18,f5.0,f5.0) endif ! Temperature @@ -1080,8 +1076,6 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si cdata_all(25,iout)=var_jb ! non linear qc if(perturb_obs) & cdata_all(26,iout)=ran01dom()*perturb_fact ! t perturbation - write(1000,1001) nread,tdiff,tvflg,pob_mb,qob,rhob,obstmp(2,1),tob,qcm,usage -1001 format(i12,f8.3,f5.0,5e25.18,f5.0,f5.0) endif ! Specific humidity if(lqob) then @@ -1117,9 +1111,6 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si cdata_all(26,iout)=var_jb ! non linear qc if(perturb_obs) & cdata_all(27,iout)=ran01dom()*perturb_fact ! q perturbation -! write(2000,1002)nread,tdiff,tvflg,pob_mb,tob,qob,rhob,qoe,obserr*one_tenth,qcm,usage -1002 format(i12,f8.3,f5.0,6e20.12,i5,1x,f5.0) - endif ! Winds --- surface wind speed diff --git a/src/gsi/setupbend.f90 b/src/gsi/setupbend.f90 index 2fe97bc072..3ff8a138cc 100644 --- a/src/gsi/setupbend.f90 +++ b/src/gsi/setupbend.f90 @@ -140,7 +140,7 @@ subroutine setupbend(obsLL,odiagLL, & use gsi_4dvar, only: nobs_bins,hr_obsbin use guess_grids, only: ges_lnprsi,hrdifsig,geop_hgti,nfldsig - use guess_grids, only: nsig_ext,gpstop + use guess_grids, only: nsig_ext,gpstop,commgpstop,commgpserrinf use gridmod, only: nsig use gridmod, only: get_ij,latlon11 use constants, only: fv,n_a,n_b,n_c,deg2rad,tiny_r_kind,r0_01,r18,r61,r63,r10000 @@ -257,6 +257,7 @@ subroutine setupbend(obsLL,odiagLL, & real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q type(obsLList),pointer,dimension(:):: gpshead + logical:: commdat gpshead => obsLL(:) save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf @@ -339,33 +340,33 @@ subroutine setupbend(obsLL,odiagLL, & nreal=mreal if (lobsdiagsave) nreal=nreal+4*miter+1 if (save_jacobian) then - nnz = nsig * 3 ! number of non-zero elements in dH(x)/dx profile - nind = 3 ! number of dense subarrays - call new(dhx_dx, nnz, nind) - nreal = nreal + size(dhx_dx) + nnz = nsig * 3 ! number of non-zero elements in dH(x)/dx profile + nind = 3 ! number of dense subarrays + call new(dhx_dx, nnz, nind) + nreal = nreal + size(dhx_dx) ! jacobian sparse array indices are the same for all obs and can be filled ! in once here: - t_ind = getindex(svars3d, 'tv') - if (t_ind < 0) then - print *, 'Error: no variable tv in state vector. Exiting.' - call stop2(1300) - endif - q_ind = getindex(svars3d, 'q') - if (q_ind < 0) then - print *, 'Error: no variable q in state vector. Exiting.' - call stop2(1300) - endif - p_ind = getindex(svars3d, 'prse') - if (p_ind < 0) then - print *, 'Error: no variable prse in state vector. Exiting.' - call stop2(1300) - endif - dhx_dx%st_ind(1) = sum(levels(1:t_ind-1)) + 1 - dhx_dx%end_ind(1) = sum(levels(1:t_ind-1)) + nsig - dhx_dx%st_ind(2) = sum(levels(1:q_ind-1)) + 1 - dhx_dx%end_ind(2) = sum(levels(1:q_ind-1)) + nsig - dhx_dx%st_ind(3) = sum(levels(1:p_ind-1)) + 1 - dhx_dx%end_ind(3) = sum(levels(1:p_ind-1)) + nsig + t_ind = getindex(svars3d, 'tv') + if (t_ind < 0) then + print *, 'Error: no variable tv in state vector. Exiting.' + call stop2(1300) + endif + q_ind = getindex(svars3d, 'q') + if (q_ind < 0) then + print *, 'Error: no variable q in state vector. Exiting.' + call stop2(1300) + endif + p_ind = getindex(svars3d, 'prse') + if (p_ind < 0) then + print *, 'Error: no variable prse in state vector. Exiting.' + call stop2(1300) + endif + dhx_dx%st_ind(1) = sum(levels(1:t_ind-1)) + 1 + dhx_dx%end_ind(1) = sum(levels(1:t_ind-1)) + nsig + dhx_dx%st_ind(2) = sum(levels(1:q_ind-1)) + 1 + dhx_dx%end_ind(2) = sum(levels(1:q_ind-1)) + nsig + dhx_dx%st_ind(3) = sum(levels(1:p_ind-1)) + 1 + dhx_dx%end_ind(3) = sum(levels(1:p_ind-1)) + nsig endif if(init_pass) call gpsrhs_alloc(is,'bend',nobs,nsig,nreal,grids_dim,nsig_ext) call gpsrhs_aliases(is) @@ -528,8 +529,8 @@ subroutine setupbend(obsLL,odiagLL, & bot_layer_SR=top_layer_SR endif endif - end do - endif + end do + endif ! locate observation in model vertical grid hob=tpdpres(i) @@ -579,290 +580,294 @@ subroutine setupbend(obsLL,odiagLL, & if(ratio_errors(i) > tiny_r_kind) then ! obs inside model grid - if (alt <= six) then - if (top_layer_SR >= 1) then ! SR exists for at least one layer. Check if obs is inside - if (tpdpres(i)==ref_rad(top_layer_SR+1)) then !obs inside model SR layer - qcfail(i)=.true. - elseif (tpdpres(i) < ref_rad(top_layer_SR+1)) then !obs below model close-to-SR layer - qcfail(i)=.true. - elseif (tpdpres(i) >= ref_rad(top_layer_SR+1) .and. tpdpres(i) <= ref_rad(top_layer_SR+2)) then !source too close - qcfail(i)=.true. - else !above - endif - endif - -! check for SR in obs, will be updated in genstats. - if ( data(igps,i) >= 0.03 .and. qc_layer_SR) then - kprof = data(iprof,i) - toss_gps_sub(kprof) = max (toss_gps_sub(kprof),data(igps,i)) - endif - endif - -! get pressure (in mb), temperature and moisture at obs location - call tintrp31(ges_lnprsi(:,:,1:nsig,:),dpressure,dlat,dlon,hob,& + if (alt <= six) then + if (top_layer_SR >= 1) then ! SR exists for at least one layer. Check if obs is inside + if (tpdpres(i)==ref_rad(top_layer_SR+1)) then !obs inside model SR layer + qcfail(i)=.true. + elseif (tpdpres(i) < ref_rad(top_layer_SR+1)) then !obs below model close-to-SR layer + qcfail(i)=.true. + elseif (tpdpres(i) >= ref_rad(top_layer_SR+1) .and. tpdpres(i) <= ref_rad(top_layer_SR+2)) then !source too close + qcfail(i)=.true. + endif + endif + +! check for SR in obs, will be updated in genstats. + if ( data(igps,i) >= 0.03_r_kind .and. qc_layer_SR) then + kprof = data(iprof,i) + toss_gps_sub(kprof) = max (toss_gps_sub(kprof),data(igps,i)) + endif + endif + +! get pressure (in mb), temperature and moisture at obs location + call tintrp31(ges_lnprsi(:,:,1:nsig,:),dpressure,dlat,dlon,hob,& dtime,hrdifsig,mype,nfldsig) - ihob=hob - k1=min(max(1,ihob),nsig) - k2=max(1,min(ihob+1,nsig)) - delz=hob-float(k1) - delz=max(zero,min(delz,one)) - trefges=tges_o(k1,i)*(one-delz)+tges_o(k2,i)*delz - qrefges=qges_o(k1)*(one-delz)+qges_o(k2)*delz !Lidia - - rdiagbuf( 6,i) = ten*exp(dpressure) ! pressure at obs location (hPa) if monotone grid - rdiagbuf(18,i) = trefges ! temperature at obs location (Kelvin) if monotone grid - rdiagbuf(21,i) = qrefges ! specific humidity at obs location (kg/kg) if monotone grid - - if (.not. qcfail(i)) then ! not SR - -! Modify error to account for representativeness error. - repe_gps=one - -! UKMET-type processing - if((data(isatid,i)==41) .or.(data(isatid,i)==722).or. & - (data(isatid,i)==723).or.(data(isatid,i)==4) .or. & - (data(isatid,i)==42) .or.(data(isatid,i)==3) .or. & - (data(isatid,i)==821).or.(data(isatid,i)==421).or. & - (data(isatid,i)==440).or.(data(isatid,i)==43) .or. & - (data(isatid,i)==5)) then + ihob=hob + k1=min(max(1,ihob),nsig) + k2=max(1,min(ihob+1,nsig)) + delz=hob-float(k1) + delz=max(zero,min(delz,one)) + trefges=tges_o(k1,i)*(one-delz)+tges_o(k2,i)*delz + qrefges=qges_o(k1)*(one-delz)+qges_o(k2)*delz !Lidia + + rdiagbuf( 6,i) = ten*exp(dpressure) ! pressure at obs location (hPa) if monotone grid + rdiagbuf(18,i) = trefges ! temperature at obs location (Kelvin) if monotone grid + rdiagbuf(21,i) = qrefges ! specific humidity at obs location (kg/kg) if monotone grid + + if (.not. qcfail(i)) then ! not SR + +! Modify error to account for representativeness error. + repe_gps=one + +! UKMET-type processing + if((data(isatid,i)==41) .or.(data(isatid,i)==722).or. & + (data(isatid,i)==723).or.(data(isatid,i)==4) .or. & + (data(isatid,i)==42) .or.(data(isatid,i)==3) .or. & + (data(isatid,i)==821).or.(data(isatid,i)==421).or. & + (data(isatid,i)==440).or.(data(isatid,i)==43) .or. & + (data(isatid,i)==5)) then + if((data(ilate,i)> r40).or.(data(ilate,i)< -r40)) then + if(alt>r12) then + repe_gps=0.19032_r_kind+0.287535_r_kind*alt-0.00260813_r_kind*alt**2 + else + repe_gps=-3.20978_r_kind+1.26964_r_kind*alt-0.0622538_r_kind*alt**2 + endif + else + if(alt>r18) then + repe_gps=-1.87788_r_kind+0.354718_r_kind*alt-0.00313189_r_kind*alt**2 + else + repe_gps=-2.41024_r_kind+0.806594_r_kind*alt-0.027257_r_kind*alt**2 + endif + endif + else +! CDAAC-type processing + if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then + if ((data(ilate,i)> r40).or.(data(ilate,i)< -r40)) then + if (alt <= 8.0_r_kind) then + repe_gps=-1.0304261_r_kind+0.3203316_r_kind*alt+0.0141337_r_kind*alt**2 + elseif (alt > 8.0_r_kind.and.alt <= r12) then + repe_gps=2.1750271_r_kind+0.0431177_r_kind*alt-0.0008567_r_kind*alt**2 + else + repe_gps=-0.3447429_r_kind+0.2829981_r_kind*alt-0.0028545_r_kind*alt**2 + endif + else + if (alt <= 4.0_r_kind) then + repe_gps=0.7285212_r_kind-1.1138755_r_kind*alt+0.2311123_r_kind*alt**2 + elseif (alt <= r18.and.alt > 4.0_r_kind) then + repe_gps=-3.3878629_r_kind+0.8691249_r_kind*alt-0.0297196_r_kind*alt**2 + else + repe_gps=-2.3875749_r_kind+0.3667211_r_kind*alt-0.0037542_r_kind*alt**2 + endif + endif + else if((data(ilate,i)> r40).or.(data(ilate,i)< -r40)) then if(alt>r12) then - repe_gps=0.19032_r_kind+0.287535_r_kind*alt-0.00260813_r_kind*alt**2 + repe_gps=-0.685627_r_kind+0.377174_r_kind*alt-0.00421934_r_kind*alt**2 else - repe_gps=-3.20978_r_kind+1.26964_r_kind*alt-0.0622538_r_kind*alt**2 + repe_gps=-3.27737_r_kind+1.20003_r_kind*alt-0.0558024_r_kind*alt**2 endif else if(alt>r18) then - repe_gps=-1.87788_r_kind+0.354718_r_kind*alt-0.00313189_r_kind*alt**2 + repe_gps=-2.73867_r_kind+0.447663_r_kind*alt-0.00475603_r_kind*alt**2 else - repe_gps=-2.41024_r_kind+0.806594_r_kind*alt-0.027257_r_kind*alt**2 + repe_gps=-3.45303_r_kind+0.908216_r_kind*alt-0.0293331_r_kind*alt**2 endif endif - else -! CDAAC-type processing - if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then - if ((data(ilate,i)> r40).or.(data(ilate,i)< -r40)) then - if (alt <= 8.0_r_kind) then - repe_gps=-1.0304261_r_kind+0.3203316_r_kind*alt+0.0141337_r_kind*alt**2 - elseif (alt > 8.0_r_kind.and.alt <= r12) then - repe_gps=2.1750271_r_kind+0.0431177_r_kind*alt-0.0008567_r_kind*alt**2 - else - repe_gps=-0.3447429_r_kind+0.2829981_r_kind*alt-0.0028545_r_kind*alt**2 - endif - else - if (alt <= 4.0_r_kind) then - repe_gps=0.7285212_r_kind-1.1138755_r_kind*alt+0.2311123_r_kind*alt**2 - elseif (alt <= r18.and.alt > 4.0_r_kind) then - repe_gps=-3.3878629_r_kind+0.8691249_r_kind*alt-0.0297196_r_kind*alt**2 - else - repe_gps=-2.3875749_r_kind+0.3667211_r_kind*alt-0.0037542_r_kind*alt**2 - endif - endif - else - if((data(ilate,i)> r40).or.(data(ilate,i)< -r40)) then - if(alt>r12) then - repe_gps=-0.685627_r_kind+0.377174_r_kind*alt-0.00421934_r_kind*alt**2 - else - repe_gps=-3.27737_r_kind+1.20003_r_kind*alt-0.0558024_r_kind*alt**2 - endif - else - if(alt>r18) then - repe_gps=-2.73867_r_kind+0.447663_r_kind*alt-0.00475603_r_kind*alt**2 - else - repe_gps=-3.45303_r_kind+0.908216_r_kind*alt-0.0293331_r_kind*alt**2 - endif - endif - endif - endif - repe_gps=exp(repe_gps) ! one/modified error in (rad-1*1E3) - repe_gps= r1em3*(one/abs(repe_gps)) ! modified error in rad - ratio_errors(i) = data(ier,i)/abs(repe_gps) - - error(i)=one/data(ier,i) ! one/original error - data(ier,i)=one/data(ier,i) ! one/original error - error_adjst(i)= ratio_errors(i)* data(ier,i) !one/adjusted error - -! Extending atmosphere above interface level nsig - d_ref_rad=ref_rad(nsig)-ref_rad(nsig-1) - do k=1,nsig_ext - ref_rad(nsig+k)=ref_rad(nsig)+ k*d_ref_rad ! extended x_i - nrefges(nsig+k,i)=nrefges(nsig+k-1,i)**2/nrefges(nsig+k-2,i) ! exended N_i - end do - -! Lagrange coefficients - ref_rad(0)=ref_rad(3) - ref_rad(nsig_up+1)=ref_rad(nsig_up-2) - do k=1,nsig_up - call setq(q_w(:,k),ref_rad(k-1:k+1),3) - enddo - - muse(i)=nint(data(iuse,i)) <= jiter - -! Get refractivity index-radius and [d(ln(n))/dx] in new grid. - intloop: do j=1,grids_dim - ref_rad_s(j)=sqrt(grid_s(j)*grid_s(j)+tpdpres(i)*tpdpres(i)) !x_j - xj(j,i)=ref_rad_s(j) - hob_s=ref_rad_s(j) - call grdcrd1(hob_s,ref_rad(1),nsig_up,1) - dbend_loc(j,i)=hob_s !location of x_j with respect to extended x_i - - - if (hob_s < rsig_up) then !obs inside the new grid - ihob=hob_s + endif -! Compute refractivity and derivative at target points -! using Lagrange interpolators - call slagdw(ref_rad(ihob-1:ihob+2),ref_rad_s(j),& - q_w(:,ihob),q_w(:,ihob+1),& - w4,dw4,4) - if(ihob==1) then - w4(4)=w4(4)+w4(1); w4(1:3)=w4(2:4);w4(4)=zero - dw4(4)=dw4(4)+dw4(1);dw4(1:3)=dw4(2:4);dw4(4)=zero - ihob=ihob+1 - endif - if(ihob==nsig_up-1) then - w4(1)=w4(1)+w4(4); w4(2:4)=w4(1:3);w4(1)=zero - dw4(1)=dw4(1)+dw4(4); dw4(2:4)=dw4(1:3);dw4(1)=zero - ihob=ihob-1 - endif - ddnj(j)=dot_product(dw4,nrefges(ihob-1:ihob+2,i))!derivative (dN/dx)_j - ddnj(j)=max(zero,abs(ddnj(j))) - else - obs_check=.true. - if (obs_check) then ! only once per obs here - nobs_out=nobs_out+1 - d_ref_rad=ref_rad(nsig)-ref_rad(nsig-1) - do kk=1,20 - ref_rad_out(nsig_up+kk)=ref_rad(nsig_up)+ kk*d_ref_rad ! extended x_i - end do - do kk=1,nsig_up - ref_rad_out(kk)=ref_rad(kk) - enddo - endif - hob_s=ref_rad_s(j) - call grdcrd1(hob_s,ref_rad_out,nsig_up+20,1) - hob_s_top=max(hob_s,hob_s_top) - endif !obs in new grid - end do intloop - - if (obs_check) then ! reject observation - qcfail(i)=.true. - data(ier,i) = zero - ratio_errors(i) = zero - muse(i)=.false. - cycle loopoverobs1 - endif - -! bending angle (radians) - dbend=ds*ddnj(1)/ref_rad_s(1) - do j=2,grids_dim - ddbend=ds*ddnj(j)/ref_rad_s(j) - dbend=dbend+two*ddbend - end do - dbend=r1em6*tpdpres(i)*dbend - -! Accumulate diagnostic information - rdiagbuf( 5,i) = (data(igps,i)-dbend)/data(igps,i) ! incremental bending angle (x100 %) - - data(igps,i)=data(igps,i)-dbend !innovation vector + repe_gps=exp(repe_gps) ! one/modified error in (rad-1*1E3) + repe_gps= r1em3*(one/abs(repe_gps)) ! modified error in rad + commdat=.false. + if (data(isatid,i)>=265 .and. data(isatid,i)<=269) then + commdat=.true. + repe_gps=commgpserrinf*repe_gps ! Inflate error for commercial data + endif + ratio_errors(i) = data(ier,i)/abs(repe_gps) + + error(i)=one/data(ier,i) ! one/original error + data(ier,i)=one/data(ier,i) ! one/original error + error_adjst(i)= ratio_errors(i)* data(ier,i) !one/adjusted error + +! Extending atmosphere above interface level nsig + d_ref_rad=ref_rad(nsig)-ref_rad(nsig-1) + do k=1,nsig_ext + ref_rad(nsig+k)=ref_rad(nsig)+ k*d_ref_rad ! extended x_i + nrefges(nsig+k,i)=nrefges(nsig+k-1,i)**2/nrefges(nsig+k-2,i) ! exended N_i + end do + +! Lagrange coefficients + ref_rad(0)=ref_rad(3) + ref_rad(nsig_up+1)=ref_rad(nsig_up-2) + do k=1,nsig_up + call setq(q_w(:,k),ref_rad(k-1:k+1),3) + enddo + + muse(i)=nint(data(iuse,i)) <= jiter + +! Get refractivity index-radius and [d(ln(n))/dx] in new grid. + intloop: do j=1,grids_dim + ref_rad_s(j)=sqrt(grid_s(j)*grid_s(j)+tpdpres(i)*tpdpres(i)) !x_j + xj(j,i)=ref_rad_s(j) + hob_s=ref_rad_s(j) + call grdcrd1(hob_s,ref_rad(1),nsig_up,1) + dbend_loc(j,i)=hob_s !location of x_j with respect to extended x_i + - if (alt <= gpstop) then ! go into qc checks - cgrossuse=cgross(ikx) - cermaxuse=cermax(ikx) - cerminuse=cermin(ikx) - if (alt > five) then - cgrossuse=cgrossuse*r400 - cermaxuse=cermaxuse*r400 - cerminuse=cerminuse*r100 + if (hob_s < rsig_up) then !obs inside the new grid + ihob=hob_s + +! Compute refractivity and derivative at target points +! using Lagrange interpolators + call slagdw(ref_rad(ihob-1:ihob+2),ref_rad_s(j),& + q_w(:,ihob),q_w(:,ihob+1),& + w4,dw4,4) + if(ihob==1) then + w4(4)=w4(4)+w4(1); w4(1:3)=w4(2:4);w4(4)=zero + dw4(4)=dw4(4)+dw4(1);dw4(1:3)=dw4(2:4);dw4(4)=zero + ihob=ihob+1 endif -! Gross error check - obserror = one/max(ratio_errors(i)*data(ier,i),tiny_r_kind) - obserrlm = max(cerminuse,min(cermaxuse,obserror)) - residual = abs(data(igps,i)) - ratio = residual/obserrlm - - if (ratio > cgrossuse) then - if (luse(i)) then - awork(4) = awork(4)+one - endif - qcfail_gross(i)=one - data(ier,i) = zero - ratio_errors(i) = zero - muse(i)=.false. - else -! Statistics QC check if obs passed gross error check - cutoff=zero - if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then - cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*one/two - else - cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*two/three - end if - cutoff2=1.5_r_kind+one*cos(data(ilate,i)*deg2rad) - if(trefges<=r240) then - cutoff3=two - else - cutoff3=0.005_r_kind*trefges**2-2.3_r_kind*trefges+266_r_kind - endif - if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then - cutoff3=cutoff3*one/two - else - cutoff3=cutoff3*two/three - end if - if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then - cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*one/two - else - cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*two/three - end if - cutoff12=((36_r_kind-alt)/two)*cutoff2+& - ((alt-34_r_kind)/two)*cutoff1 - cutoff23=((eleven-alt)/two)*cutoff3+& - ((alt-nine)/two)*cutoff2 - cutoff34=((six-alt)/two)*cutoff4+& - ((alt-four)/two)*cutoff3 - if(alt>36_r_kind) cutoff=cutoff1 - if((alt<=36_r_kind).and.(alt>34_r_kind)) cutoff=cutoff12 - if((alt<=34_r_kind).and.(alt>eleven)) cutoff=cutoff2 - if((alt<=eleven).and.(alt>nine)) cutoff=cutoff23 - if((alt<=nine).and.(alt>six)) cutoff=cutoff3 - if((alt<=six).and.(alt>four)) cutoff=cutoff34 - if(alt<=four) cutoff=cutoff4 - - if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then - cutoff=two*cutoff*r0_01 - else - cutoff=three*cutoff*r0_01 - end if - - if(abs(rdiagbuf(5,i)) > cutoff) then - qcfail(i)=.true. - data(ier,i) = zero - ratio_errors(i) = zero - muse(i) = .false. - end if - end if ! gross qc check - end if ! qc checks (only below 50km) - -! Remove obs above 50 km - if(alt > gpstop) then - data(ier,i) = zero - ratio_errors(i) = zero - qcfail_high(i)=one - muse(i)=.false. - endif + if(ihob==nsig_up-1) then + w4(1)=w4(1)+w4(4); w4(2:4)=w4(1:3);w4(1)=zero + dw4(1)=dw4(1)+dw4(4); dw4(2:4)=dw4(1:3);dw4(1)=zero + ihob=ihob-1 + endif + ddnj(j)=dot_product(dw4,nrefges(ihob-1:ihob+2,i))!derivative (dN/dx)_j + ddnj(j)=max(zero,abs(ddnj(j))) + else + obs_check=.true. + if (obs_check) then ! only once per obs here + nobs_out=nobs_out+1 + d_ref_rad=ref_rad(nsig)-ref_rad(nsig-1) + do kk=1,20 + ref_rad_out(nsig_up+kk)=ref_rad(nsig_up)+ kk*d_ref_rad ! extended x_i + end do + do kk=1,nsig_up + ref_rad_out(kk)=ref_rad(kk) + enddo + endif + hob_s=ref_rad_s(j) + call grdcrd1(hob_s,ref_rad_out,nsig_up+20,1) + hob_s_top=max(hob_s,hob_s_top) + endif !obs in new grid + end do intloop + + if (obs_check) then ! reject observation + qcfail(i)=.true. + data(ier,i) = zero + ratio_errors(i) = zero + muse(i)=.false. + cycle loopoverobs1 + endif + +! bending angle (radians) + dbend=ds*ddnj(1)/ref_rad_s(1) + do j=2,grids_dim + ddbend=ds*ddnj(j)/ref_rad_s(j) + dbend=dbend+two*ddbend + end do + dbend=r1em6*tpdpres(i)*dbend + +! Accumulate diagnostic information + rdiagbuf( 5,i) = (data(igps,i)-dbend)/data(igps,i) ! incremental bending angle (x100 %) + + data(igps,i)=data(igps,i)-dbend !innovation vector + if (alt <= gpstop) then ! go into qc checks + if ((alt <= commgpstop) .or. (.not.commdat)) then + cgrossuse=cgross(ikx) + cermaxuse=cermax(ikx) + cerminuse=cermin(ikx) + if (alt > five) then + cgrossuse=cgrossuse*r400 + cermaxuse=cermaxuse*r400 + cerminuse=cerminuse*r100 + endif +! Gross error check + obserror = one/max(ratio_errors(i)*data(ier,i),tiny_r_kind) + obserrlm = max(cerminuse,min(cermaxuse,obserror)) + residual = abs(data(igps,i)) + ratio = residual/obserrlm + + if (ratio > cgrossuse) then + if (luse(i)) then + awork(4) = awork(4)+one + endif + qcfail_gross(i)=one + data(ier,i) = zero + ratio_errors(i) = zero + muse(i)=.false. + else +! Statistics QC check if obs passed gross error check + cutoff=zero + if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then + cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*one/two + else + cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*two/three + end if + cutoff2=1.5_r_kind+one*cos(data(ilate,i)*deg2rad) + if(trefges<=r240) then + cutoff3=two + else + cutoff3=0.005_r_kind*trefges**2-2.3_r_kind*trefges+266_r_kind + endif + if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then + cutoff3=cutoff3*one/two + else + cutoff3=cutoff3*two/three + end if + if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then + cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*one/two + else + cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*two/three + end if + cutoff12=((36_r_kind-alt)/two)*cutoff2+& + ((alt-34_r_kind)/two)*cutoff1 + cutoff23=((eleven-alt)/two)*cutoff3+& + ((alt-nine)/two)*cutoff2 + cutoff34=((six-alt)/two)*cutoff4+& + ((alt-four)/two)*cutoff3 + if(alt>36_r_kind) cutoff=cutoff1 + if((alt<=36_r_kind).and.(alt>34_r_kind)) cutoff=cutoff12 + if((alt<=34_r_kind).and.(alt>eleven)) cutoff=cutoff2 + if((alt<=eleven).and.(alt>nine)) cutoff=cutoff23 + if((alt<=nine).and.(alt>six)) cutoff=cutoff3 + if((alt<=six).and.(alt>four)) cutoff=cutoff34 + if(alt<=four) cutoff=cutoff4 + + if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then + cutoff=two*cutoff*r0_01 + else + cutoff=three*cutoff*r0_01 + end if + + if(abs(rdiagbuf(5,i)) > cutoff) then + qcfail(i)=.true. + data(ier,i) = zero + ratio_errors(i) = zero + muse(i) = .false. + end if + end if !gross qc check + end if ! commdat < commgpstop + end if ! qc checks (only below 50km) +! Remove obs above 50 km + if((alt > gpstop) .or. (commdat .and. (alt > commgpstop))) then + data(ier,i) = zero + ratio_errors(i) = zero + qcfail_high(i)=one + muse(i)=.false. + endif + +! Remove MetOP/GRAS data below 8 km + if( (alt <= eight) .and. & + ((data(isatid,i)==4).or.(data(isatid,i)==3).or.(data(isatid,i)==5))) then + qcfail(i)=.true. + data(ier,i) = zero + ratio_errors(i) = zero + muse(i)=.false. + endif -! Remove MetOP/GRAS data below 8 km - if( (alt <= eight) .and. & - ((data(isatid,i)==4).or.(data(isatid,i)==3).or.(data(isatid,i)==5))) then - qcfail(i)=.true. - data(ier,i) = zero - ratio_errors(i) = zero - muse(i)=.false. - endif - - end if ! obs above super-refraction and shadow layers + end if ! obs above super-refraction and shadow layers end if ! obs inside the vertical grid end do loopoverobs1 ! end of loop over observations @@ -1056,7 +1061,7 @@ subroutine setupbend(obsLL,odiagLL, & allocate(my_head) call gpsNode_appendto(my_head,gpshead(ibin)) - + my_head%idv = is my_head%iob = ioid(i) my_head%elat= data(ilate,i) @@ -1199,7 +1204,7 @@ subroutine setupbend(obsLL,odiagLL, & data(ier,i) = zero rdiagbuf(12,i) = -one rdiagbuf(10,i) = six - end if + end if if (save_jacobian) then @@ -1279,55 +1284,55 @@ subroutine init_vars_ varname='z' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) if (istatus==0) then - if(allocated(ges_z))then - write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' - call stop2(999) - endif - allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig)) - ges_z(:,:,1)=rank2 - do ifld=2,nfldsig - call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) - ges_z(:,:,ifld)=rank2 - enddo + if(allocated(ges_z))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig)) + ges_z(:,:,1)=rank2 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) + ges_z(:,:,ifld)=rank2 + enddo else - write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus - call stop2(999) + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) endif ! get tv ... varname='tv' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) if (istatus==0) then - if(allocated(ges_tv))then - write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' - call stop2(999) - endif - allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) - ges_tv(:,:,:,1)=rank3 - do ifld=2,nfldsig - call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) - ges_tv(:,:,:,ifld)=rank3 - enddo + if(allocated(ges_tv))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_tv(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_tv(:,:,:,ifld)=rank3 + enddo else - write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus - call stop2(999) + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) endif ! get q ... varname='q' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) if (istatus==0) then - if(allocated(ges_q))then - write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' - call stop2(999) - endif - allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) - ges_q(:,:,:,1)=rank3 - do ifld=2,nfldsig - call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) - ges_q(:,:,:,ifld)=rank3 - enddo + if(allocated(ges_q))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_q(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_q(:,:,:,ifld)=rank3 + enddo else - write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus - call stop2(999) + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) endif else write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index 7f001ea9cd..a5b2a6b285 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -2108,7 +2108,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& deallocate(diagbufchan) deallocate(sc_index) - if (rad_diagsave) then + if (rad_diagsave .and. nchanl_diag > 0) then if (netcdf_diag) call nc_diag_write if(binary_diag) call final_binary_diag_ if (lextra .and. allocated(diagbufex)) deallocate(diagbufex) diff --git a/src/gsi/stop1.f90 b/src/gsi/stop1.f90 index 3682f3bdee..39f853fc61 100644 --- a/src/gsi/stop1.f90 +++ b/src/gsi/stop1.f90 @@ -23,13 +23,20 @@ subroutine stop2(ierror_code) ! !$$$ use kinds, only: i_kind - use mpimod, only: mpi_comm_world,ierror + use mpimod, only: mpi_comm_world,ierror,mype + use gsi_io, only: verbose implicit none integer(i_kind) ierror_code - write(6,*)'****STOP2**** ABORTING EXECUTION w/code=',ierror_code - write(0,*)'****STOP2**** ABORTING EXECUTION w/code=',ierror_code + if (verbose) then + write(6,*)'****STOP2**** ABORTING EXECUTION w/code=',ierror_code + write(0,*)'****STOP2**** ABORTING EXECUTION w/code=',ierror_code + elseif (mype==0) then + write(6,*)'****STOP2**** ABORTING EXECUTION w/code=',ierror_code + write(0,*)'****STOP2**** ABORTING EXECUTION w/code=',ierror_code + endif + call mpi_abort(mpi_comm_world,ierror_code,ierror) stop return diff --git a/src/gsi/write_incr.f90 b/src/gsi/write_incr.f90 index eae5c5ffe4..6563775662 100644 --- a/src/gsi/write_incr.f90 +++ b/src/gsi/write_incr.f90 @@ -253,14 +253,14 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) ! end the netCDF file definition call nccheck_incr(nf90_enddef(ncid_out)) - ! compute delz + ! compute delz (so that delz < 0 as model expects) do k=1,grd%nsig - sub_dzb(:,:,k) = ges_geopi(:,:,k+1,ibin) - ges_geopi(:,:,k,ibin) + sub_dzb(:,:,k) = ges_geopi(:,:,k,ibin) - ges_geopi(:,:,k+1,ibin) enddo call load_geop_hgt do k=1,grd%nsig - sub_dza(:,:,k) = geop_hgti(:,:,k+1,ibin) - geop_hgti(:,:,k,ibin) + sub_dza(:,:,k) = geop_hgti(:,:,k,ibin) - geop_hgti(:,:,k+1,ibin) enddo sub_dza = sub_dza - sub_dzb !sub_dza is increment diff --git a/ush/build_all_cmake.sh b/ush/build_all_cmake.sh index 02f46cf7b5..8f32064872 100755 --- a/ush/build_all_cmake.sh +++ b/ush/build_all_cmake.sh @@ -7,7 +7,17 @@ pwd=$(pwd) build_type=${1:-'PRODUCTION'} dir_root=${2:-$pwd} +mode=${3:-'EMC'} + +# If NCO build, prune directories and files before build +if [ $mode = NCO ]; then + cd $dir_root/ush + $dir_root/ush/prune_4nco_global.sh prune +fi + + +# Initialize and load modules if [[ -d /dcom && -d /hwrf ]] ; then . /usrx/local/Modules/3.2.10/init/sh target=wcoss @@ -89,6 +99,17 @@ else cmake .. fi -make -j 8 +# Build apps. Echo extra printout for NCO build +if [ $mode = NCO ]; then + make VERBOSE=1 -j 8 +else + make -j 8 +fi +rc=$? + +# If NCO build is successful, remove build directory +if [ $mode = NCO -a $rc -eq 0 ]; then + rm -rf $dir_root/build +fi exit diff --git a/ush/prune_4nco_global.sh b/ush/prune_4nco_global.sh new file mode 100755 index 0000000000..699064a613 --- /dev/null +++ b/ush/prune_4nco_global.sh @@ -0,0 +1,189 @@ +#!/bin/sh + +# This script removes or restores directories and/or files found in +# the rlist below from the current working directory. This script +# is intended to be executed when NCO implements GFS DA updates. +# NCO does not want package installations for operations to include +# non-operational code. This script removes directories and files +# not used by operations from the installation directory. +# +# Two modes are supported: prune, restore +# prune: use git rm -r to remove directories and files. +# The directories and files to remove are found +# in rlist below +# restore: use git RESET head and git checkout to restore +# removed directories and files +# + +set -ex + +mode=$1 + +# Check mode and set string +if [[ "$mode" = "prune" ]]; then + string="rm -r" +elif [[ "$mode" = "restore" ]]; then + string="reset HEAD" +else + echo " " + echo "***ERROR*** invalid mode= $mode" + echo " valid modes are prune or restore" + echo " " + exit +fi + + +# Set root directory +cd .. +topdir=$(pwd) + +echo " " +echo "Execute git $string in $topdir" +echo " " + + +# Process top level directories +cd $topdir +rlist="regression src/GSD unit-tests" +for type in $rlist; do + git $string ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git $string ${type}" + exit + fi + if [[ "$mode" = "restore" ]]; then + git checkout ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git checkout ${type}" + exit + fi + fi +done + + +# Process doc directories and files +cd $topdir/doc +rlist="EnKF_user_guide GSI_user_guide README.discover" +for type in $rlist; do + git $string ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git $string ${type}" + exit + fi + if [[ "$mode" = "restore" ]]; then + git checkout ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git checkout ${type}" + exit + fi + fi +done + + +# Process scripts directories and files +cd $topdir/scripts +rlist="exurma2p5_gsianl.sh" +for type in $rlist; do + git $string ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git $string ${type}" + exit + fi + if [[ "$mode" = "restore" ]]; then + git checkout ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git checkout ${type}" + exit + fi + fi +done + + +# Process ush directories and files +cd $topdir/ush +rlist="build.comgsi EnKF com Get_Initial_Files gfs_truncate_enkf llsub para refactor_4nco_global run_arw rungsi sub" +for type in $rlist; do + git $string ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git $string ${type}" + exit + fi + if [[ "$mode" = "restore" ]]; then + git checkout ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git checkout ${type}" + exit + fi + fi +done + + +# Process util directories and files +cd $topdir/util +rlist="Aero Analysis_Utilities Baseline Config Correlated_Obs DTC EFSOI FOV GEN_BE_V2.0 GMI_BUFR MODIS_AOD Misc NCEP NMC_Bkerror README Radar_Monitor Radiance_bias_correction_Utilities Radiance_Utilities Single_Observation bufr_tools global_angupdate gsienvreport.sh ndate python_utilities radar_process zero_biascoeff" +for type in $rlist; do + git $string ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git $string ${type}" + exit + fi + if [[ "$mode" = "restore" ]]; then + git checkout ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git checkout ${type}" + exit + fi + fi +done + + +# Process util/EnKF directories and files +cd $topdir/util/EnKF +rlist="arw python_utilities" +for type in $rlist; do + git $string ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git $string ${type}" + exit + fi + if [[ "$mode" = "restore" ]]; then + git checkout ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git checkout ${type}" + exit + fi + fi +done + + +# Process util/EnKF/gfs/src directories and files +cd $topdir/util/EnKF/gfs/src +rlist="adjustps misc preproc gribmean recenterncio_hybgain recenternemsiop_hybgain" +for type in $rlist; do + git $string ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git $string ${type}" + exit + fi + if [[ "$mode" = "restore" ]]; then + git checkout ${type}* + rc=$? + if [[ $rc -ne 0 ]]; then + echo "***ERROR* git checkout ${type}" + exit + fi + fi +done diff --git a/util/Minimization_Monitor/nwprod/gdas.v1.0.0/scripts/exgdas_atmos_vminmon.sh b/util/Minimization_Monitor/nwprod/gdas.v1.0.0/scripts/exgdas_atmos_vminmon.sh index af3e338151..1fa612cbfe 100755 --- a/util/Minimization_Monitor/nwprod/gdas.v1.0.0/scripts/exgdas_atmos_vminmon.sh +++ b/util/Minimization_Monitor/nwprod/gdas.v1.0.0/scripts/exgdas_atmos_vminmon.sh @@ -59,7 +59,6 @@ export mm_costfile=${costfile:-${M_FIXgdas}/gdas_minmon_cost.txt} export MINMON_SUFFIX=${MINMON_SUFFIX:-GDAS} export PDATE=${PDY}${cyc} export NCP=${NCP:-/bin/cp} -export NDATE=${NDATE:-/nwprod/util/exec/ndate} export pgm=exgdas_vrfminmon.sh if [[ ! -d ${DATA} ]]; then diff --git a/util/Minimization_Monitor/nwprod/gfs.v1.0.0/scripts/exgfs_atmos_vminmon.sh b/util/Minimization_Monitor/nwprod/gfs.v1.0.0/scripts/exgfs_atmos_vminmon.sh index 98e33c708d..eb0eac23c5 100755 --- a/util/Minimization_Monitor/nwprod/gfs.v1.0.0/scripts/exgfs_atmos_vminmon.sh +++ b/util/Minimization_Monitor/nwprod/gfs.v1.0.0/scripts/exgfs_atmos_vminmon.sh @@ -58,7 +58,6 @@ export mm_costfile=${costfile:-${M_FIXgfs}/gfs_minmon_cost.txt} export MINMON_SUFFIX=${MINMON_SUFFIX:-GFS} export PDATE=${PDY}${cyc} export NCP=${NCP:-/bin/cp} -export NDATE=${NDATE:-/nwprod/util/exec/ndate} export pgm=exgfs_vrfminmon.sh diff --git a/util/Radiance_Monitor/nwprod/gdas_radmon/scripts/exgdas_atmos_verfrad.sh b/util/Radiance_Monitor/nwprod/gdas_radmon/scripts/exgdas_atmos_verfrad.sh index 06244e2c04..263c3b59a6 100755 --- a/util/Radiance_Monitor/nwprod/gdas_radmon/scripts/exgdas_atmos_verfrad.sh +++ b/util/Radiance_Monitor/nwprod/gdas_radmon/scripts/exgdas_atmos_verfrad.sh @@ -57,7 +57,6 @@ export USE_MAIL=${USE_MAIL:-0} export MAIL_TO=${MAIL_TO:-" "} export MAIL_CC=${MAIL_CC:-" "} export NCP=${NCP:-/bin/cp} -export NDATE=${NDATE:-/nwprod/util/exec/ndate} ########################################################################### # ensure TANK dir exists, verify radstat and biascr are available @@ -208,6 +207,15 @@ elif [[ $rc_time -ne 0 ]]; then err=$rc_time fi +##################################################################### +# Restrict select sensors and satellites +export CHGRP_CMD=${CHGRP_CMD:-"chgrp ${group_name:-rstprod}"} +rlist="saphir" +for rtype in $rlist; do + ${CHGRP_CMD} $TANKverf_rad/*${rtype}* +done + + if [[ "$VERBOSE" = "YES" ]]; then echo "end exgdas_vrfyrad.sh, exit value = ${err}" fi