diff --git a/.gitmodules b/.gitmodules index ee727b809a..f4bfd21ee7 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,7 @@ +[submodule "fix"] + path = fix + url = gerrit:GSI-fix + [submodule "libsrc"] path = libsrc url = gerrit:GSI-libsrc diff --git a/CMakeLists.txt b/CMakeLists.txt index 38310a7f2a..2603bc76d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -61,7 +61,7 @@ project(GSI) option(BUILD_WRF "Build the Enkf with WRF module " OFF) endif() option(BUILD_UTIL "Build the Enkf utilities " OFF) - option(BUILD_EFSOI "Build the EFSOI utility " OFF) + option(BUILD_UTIL_COM "Build community utilities " OFF) option(BUILD_ENKF_PREPROCESS_ARW "Build enkf preprocess for ARW " OFF) option(BUILD_COV_CALC "Build the Desroziers utility" OFF) @@ -273,10 +273,7 @@ project(GSI) add_subdirectory(util/Radiance_Monitor/nwprod/radmon_shared.v3.0.0/sorc/verf_radbcoef.fd) add_subdirectory(util/Radiance_Monitor/nwprod/radmon_shared.v3.0.0/sorc/verf_radbcor.fd) add_subdirectory(util/Radiance_Monitor/nwprod/radmon_shared.v3.0.0/sorc/verf_radtime.fd) - endif(BUILD_UTIL) - if(BUILD_EFSOI) - add_subdirectory(util/EFSOI_Utilities/src) - endif(BUILD_EFSOI) + endif(BUILD_UTIL) find_package( NDATE ) if( NOT NDATE ) add_subdirectory(util/ndate) diff --git a/jobs/JGDAS_EFSOI b/jobs/JGDAS_EFSOI deleted file mode 100755 index dbd8ab86f8..0000000000 --- a/jobs/JGDAS_EFSOI +++ /dev/null @@ -1,142 +0,0 @@ -#!/bin/ksh -set -x - -export RUN_ENVIR=${RUN_ENVIR:-"nco"} -export PS4='$SECONDS + ' -date - - -############################# -# Source relevant config files -############################# -export EXPDIR=${EXPDIR:-$HOMEgfs/parm/config} - -configs="base anal efsoi" -config_path=${EXPDIR:-$NWROOT/gfs.${gfs_ver}/parm/config} -for config in $configs; do - . $config_path/config.$config - status=$? - [[ $status -ne 0 ]] && exit $status -done - - -########################################## -# Source machine runtime environment -########################################## -. $HOMEgfs/env/${machine}.env efsoi - - -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# Obtain unique process id (pid) and make temp directory -############################################## -export pid=${pid:-$$} -export outid=${outid:-"LL$job"} -export DATA=${DATA:-${DATAROOT}/${jobid:?}} -mkdir -p $DATA -cd $DATA - - -############################################## -# Run setpdy and initialize PDY variables -############################################## -export cycle="t${cyc}z" -setpdy.sh -. ./PDY - - -############################################## -# Determine Job Output Name on System -############################################## -export pgmout="OUTPUT.${pid}" -export pgmerr=errfile - - -############################################## -# Set variables used in the exglobal script -############################################## -export CDATE=${CDATE:-${PDY}${cyc}} -export CDUMP=${CDUMP:-${RUN:-"gdas"}} -export COMPONENT=${COMPONENT:-atmos} -if [ $RUN_ENVIR = "nco" ]; then - export ROTDIR=${COMROOT:?}/$NET/$envir -fi - - -############################################## -# Begin JOB SPECIFIC work -############################################## - -GDATE=$($NDATE -$assim_freq $CDATE) -gPDY=$(echo $GDATE | cut -c1-8) -export gcyc=$(echo $GDATE | cut -c9-10) # needed in namelist - -VDATE=$($NDATE 24 $CDATE) -vPDY=$(echo $VDATE | cut -c1-8) -vcyc=$(echo $VDATE | cut -c9-10) - -export APREFIX="${CDUMP}.t${cyc}z." -export GPREFIX="gdas.t${gcyc}z." -export ASUFFIX=${ASUFFIX:-$SUFFIX} -export GSUFFIX=${GSUFFIX:-$SUFFIX} - - - -# COMIN_GES_ENS and COMOUT_ANL_ENS are used in script -export COMIN_GES_ENS="$ROTDIR/efsoigdas.$gPDY/$gcyc/$COMPONENT" -export COMOUT_ANL_ENS="$ROTDIR/enkf$CDUMP.$PDY/$cyc/$COMPONENT" -export COMIN_ANL="$ROTDIR/$CDUMP.$vPDY/$vcyc/$COMPONENT" -export COMOUT_ANL_ENSFSOI="$ROTDIR/efsoi$CDUMP.$PDY/$cyc/$COMPONENT" -export OSENSE_SAVE_DIR="$ROTDIR/osense" -mkdir -p $OSENSE_SAVE_DIR -# ======================================================== - -############################################################### -# Run relevant exglobal script -env -msg="HAS BEGUN on `hostname`" -postmsg "$jlogfile" "$msg" -$LOGSCRIPT - -${EFSOIUPDSH:-$SCRgfs/exgdas_efsoi.sh} -status=$? -[[ $status -ne 0 ]] && exit $status - - - -############################################## -# Send Alerts -############################################## -if [ $SENDDBN = YES ] ; then - $DBNROOT/bin/dbn_alert MODEL ENKF1_MSC_enkfstat $job $COMOUT_ANL_ENSFSOI/${APREFIX}enkfstat -fi - - -############################################## -# End JOB SPECIFIC work -############################################## - -############################################## -# Final processing -############################################## -if [ -e "$pgmout" ] ; then - cat $pgmout -fi - - -msg="ENDED NORMALLY." -postmsg "$jlogfile" "$msg" - - -########################################## -# Remove the Temporary working directory -########################################## -cd $DATAROOT - -[[ $KEEPDATA = "NO" ]] && rm -rf $DATA - -date -exit 0 diff --git a/jobs/JGDAS_EFSOI_ECEN b/jobs/JGDAS_EFSOI_ECEN deleted file mode 100755 index 28d04b9f5c..0000000000 --- a/jobs/JGDAS_EFSOI_ECEN +++ /dev/null @@ -1,143 +0,0 @@ -#!/bin/ksh -set -x - -export RUN_ENVIR=${RUN_ENVIR:-"nco"} -export PS4='$SECONDS + ' -date - - -############################# -# Source relevant config files -############################# -export EXPDIR=${EXPDIR:-$HOMEgfs/parm/config} -configs="base ecenfsoi" -config_path=${EXPDIR:-$NWROOT/gfs.${gfs_ver}/parm/config} -for config in $configs; do - . $config_path/config.$config - status=$? - [[ $status -ne 0 ]] && exit $status -done - - -########################################## -# Source machine runtime environment -########################################## -. $HOMEgfs/env/${machine}.env ecenfsoi -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# Obtain unique process id (pid) and make temp directory -############################################## -export pid=${pid:-$$} -export outid=${outid:-"LL$job"} -export DATA=${DATA:-${DATAROOT}/${jobid:?}} -mkdir -p $DATA -cd $DATA - - -############################################## -# Run setpdy and initialize PDY variables -############################################## -export cycle="t${cyc}z" -setpdy.sh -. ./PDY - - -############################################## -# Determine Job Output Name on System -############################################## -export pgmout="OUTPUT.${pid}" -export pgmerr=errfile - - -############################################## -# Set variables used in the script -############################################## -export CDATE=${CDATE:-${PDY}${cyc}} -export CDUMP=${CDUMP:-${RUN:-"gdas"}} -export COMPONENT=${COMPONENT:-atmos} -if [ $RUN_ENVIR = "nco" ]; then - export ROTDIR=${COMROOT:?}/$NET/$envir -fi - - -############################################## -# Begin JOB SPECIFIC work -############################################## - -GDATE=$($NDATE -$assim_freq $CDATE) -gPDY=$(echo $GDATE | cut -c1-8) -gcyc=$(echo $GDATE | cut -c9-10) -GDUMP=${GDUMP:-"gdas"} - -export CASE=$CASE_ENKF - - -EUPD_CYC=$(echo ${EUPD_CYC:-"gdas"} | tr a-z A-Z) -if [ $EUPD_CYC = "GFS" ]; then - CDUMP_ENKF="gfs" -else - CDUMP_ENKF=$CDUMP -fi - -export OPREFIX="${CDUMP}.t${cyc}z." -export APREFIX="${CDUMP}.t${cyc}z." -export APREFIX_ENKF="${CDUMP_ENKF}.t${cyc}z." -export GPREFIX="${CDUMP}.t${gcyc}z." -export GSUFFIX=${GSUFFIX:-$SUFFIX} -export ASUFFIX=${ASUFFIX:-$SUFFIX} - -if [ $RUN_ENVIR = "nco" -o ${ROTDIR_DUMP:-NO} = "YES" ]; then - export COMIN_OBS=${COMIN_OBS:-$ROTDIR/$RUN.$PDY/$cyc/$COMPONENT} - export COMIN_GES_OBS=${COMIN_GES_OBS:-$ROTDIR/$GDUMP.$gPDY/$gcyc/$COMPONENT} -else - export COMIN_OBS="$DMPDIR/$CDATE/$CDUMP" - export COMIN_GES_OBS="$DMPDIR/$GDATE/$GDUMP" -fi - -# COMIN, COMIN_ENS and COMIN_GES_ENS are used in script -export COMIN="$ROTDIR/$CDUMP.$PDY/$cyc/$COMPONENT" -export COMIN_ENS="$ROTDIR/efsoi$CDUMP_ENKF.$PDY/$cyc/$COMPONENT" -export COMOUT_ENS="$ROTDIR/efsoi$CDUMP.$PDY/$cyc/$COMPONENT" -export COMIN_GES_ENS="$ROTDIR/enkf$CDUMP.$gPDY/$gcyc/$COMPONENT" - - -############################################################### -# Run relevant script -env -msg="HAS BEGUN on `hostname`" -postmsg "$jlogfile" "$msg" -$LOGSCRIPT - - -${ENKFRECENSH:-$SCRgfs/exgdas_enkf_ecen.sh} -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# End JOB SPECIFIC work -############################################## - -############################################## -# Final processing -############################################## -if [ -e "$pgmout" ] ; then - cat $pgmout -fi - - -msg="ENDED NORMALLY." -postmsg "$jlogfile" "$msg" - - -########################################## -# Remove the Temporary working directory -########################################## -cd $DATAROOT -[[ $KEEPDATA = "NO" ]] && rm -rf $DATA - -date -exit 0 diff --git a/jobs/JGDAS_EFSOI_FCST b/jobs/JGDAS_EFSOI_FCST deleted file mode 100755 index 73f11f479d..0000000000 --- a/jobs/JGDAS_EFSOI_FCST +++ /dev/null @@ -1,149 +0,0 @@ -#!/bin/ksh -set -x - -export RUN_ENVIR=${RUN_ENVIR:-"nco"} -export PS4='$SECONDS + ' -date - - -############################# -# Source relevant config files -############################# -export EXPDIR=${EXPDIR:-$HOMEgfs/parm/config} -config_path=${EXPDIR:-$NWROOT/gfs.${gfs_ver}/parm/config} -configs="base fcst efcsfsoi" -for config in $configs; do - . $config_path/config.$config - status=$? - [[ $status -ne 0 ]] && exit $status -done - - -########################################## -# Source machine runtime environment -########################################## -. $HOMEgfs/env/${machine}.env efcsfsoi -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# Obtain unique process id (pid) and make temp directory -############################################## -export pid=${pid:-$$} -export outid=${outid:-"LL$job"} -export DATA=${DATA:-${DATAROOT}/${jobid:?}} -mkdir -p $DATA -cd $DATA - - -############################################## -# Run setpdy and initialize PDY variables -############################################## -export cycle="t${cyc}z" -setpdy.sh -. ./PDY - - -############################################## -# Determine Job Output Name on System -############################################## -export pgmout="OUTPUT.${pid}" -export pgmerr=errfile - - -############################################## -# Set variables used in the script -############################################## -export CDATE=${CDATE:-${PDY}${cyc}} -export CDUMP=${CDUMP:-${RUN:-"gdas"}} -export COMPONENT=${COMPONENT:-atmos} -if [ $RUN_ENVIR = "nco" ]; then - export ROTDIR=${COMROOT:?}/$NET/$envir -fi - - -############################################## -# Begin JOB SPECIFIC work -############################################## - -export CASE=$CASE_ENKF - -# COMOUT is used in script -export COMOUT="$ROTDIR/efsoi$CDUMP.$PDY/$cyc/$COMPONENT" - - -# Forecast length for EnKF forecast -export FHMIN_ENKF=${FHMIN_EFSOI:-0} -export FHOUT_ENKF=${FHOUT_EFSOI:-6} -export FHMAX_ENKF=${FHMAX_EFSOI:-30} -export FHMIN=$FHMIN_ENKF -export FHOUT=$FHOUT_ENKF -export FHMAX=$FHMAX_ENKF - -export FORECASTSH=$HOMEgfs/scripts/exglobal_efsoi_forecast.sh - -# Get ENSBEG/ENSEND from ENSGRP and NMEM_EFCSGRP -export ENSEND=$((NMEM_EFCSGRP * ENSGRP)) -export ENSBEG=$((ENSEND - NMEM_EFCSGRP + 1)) - - -############################################################### -# Run relevant script -env -msg="HAS BEGUN on `hostname`" -postmsg "$jlogfile" "$msg" -$LOGSCRIPT - - -${ENKFFCSTSH:-$SCRgfs/exgdas_enkf_fcst.sh} -status=$? -[[ $status -ne 0 ]] && exit $status - - -# Double check the status of members in ENSGRP -EFCSGRP=$COMOUT/efcs.grp${ENSGRP} -npass=0 -if [ -f $EFCSGRP ]; then - npass=$(grep "PASS" $EFCSGRP | wc -l) -fi -echo "$npass/$NMEM_EFCSGRP members successfull in efcs.grp$ENSGRP" -if [ $npass -ne $NMEM_EFCSGRP ]; then - echo "FATAL ERROR: Failed members in group $ENSGRP, ABORT!" - cat $EFCSGRP - exit 99 -fi - - -############################################## -# Send Alerts -############################################## -if [ $SENDDBN = YES ] ; then - $DBNROOT/bin/dbn_alert MODEL ENKF1_MSC_fcsstat $job $EFCSGRP -fi - - -############################################## -# End JOB SPECIFIC work -############################################## - -############################################## -# Final processing -############################################## -if [ -e "$pgmout" ] ; then - cat $pgmout -fi - - -msg="ENDED NORMALLY." -postmsg "$jlogfile" "$msg" - - -########################################## -# Remove the Temporary working directory -########################################## -cd $DATAROOT -[[ $KEEPDATA = "NO" ]] && rm -rf $DATA - -date -exit 0 diff --git a/jobs/JGDAS_EFSOI_POST b/jobs/JGDAS_EFSOI_POST deleted file mode 100755 index f71f865b7e..0000000000 --- a/jobs/JGDAS_EFSOI_POST +++ /dev/null @@ -1,121 +0,0 @@ -#!/bin/ksh -set -x - -export RUN_ENVIR=${RUN_ENVIR:-"nco"} -export PS4='$SECONDS + ' -date - - -############################# -# Source relevant config files -############################# -export EXPDIR=${EXPDIR:-$HOMEgfs/parm/config} -configs="base eposfsoi" -config_path=${EXPDIR:-$NWROOT/gfs.${gfs_ver}/parm/config} -for config in $configs; do - . $config_path/config.$config - status=$? - [[ $status -ne 0 ]] && exit $status -done - - -########################################## -# Source machine runtime environment -########################################## -. $HOMEgfs/env/${machine}.env eposfsoi -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# Obtain unique process id (pid) and make temp directory -############################################## -export pid=${pid:-$$} -export outid=${outid:-"LL$job"} -export DATA=${DATA:-${DATAROOT}/${jobid:?}} -mkdir -p $DATA -cd $DATA - - -############################################## -# Run setpdy and initialize PDY variables -############################################## -export cycle="t${cyc}z" -setpdy.sh -. ./PDY - - -############################################## -# Determine Job Output Name on System -############################################## -export pgmout="OUTPUT.${pid}" -export pgmerr=errfile - - -############################################## -# Set variables used in the script -############################################## -export CDATE=${CDATE:-${PDY}${cyc}} -export CDUMP=${CDUMP:-${RUN:-"gdas"}} -export COMPONENT=${COMPONENT:-atmos} -if [ $RUN_ENVIR = "nco" ]; then - export ROTDIR=${COMROOT:?}/$NET/$envir -fi - - -############################################## -# Begin JOB SPECIFIC work -############################################## -export GFS_NCIO=${GFS_NCIO:-"YES"} - -export PREFIX="${CDUMP}.t${cyc}z." - -# COMIN, COMOUT are used in script -export COMIN="$ROTDIR/efsoi$CDUMP.$PDY/$cyc/$COMPONENT" -export COMOUT="$ROTDIR/efsoi$CDUMP.$PDY/$cyc/$COMPONENT" - - -export LEVS=$((LEVS-1)) - -#export FHMIN_EPOS=6 -#export FHMAX_EPOS=30 -#export FHOUT_EPOS=6 - - -############################################################### -# Run relevant script -env -msg="HAS BEGUN on `hostname`" -postmsg "$jlogfile" "$msg" -$LOGSCRIPT - - -${ENKFPOSTSH:-$SCRgfs/exgdas_enkf_post.sh} -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# End JOB SPECIFIC work -############################################## - -############################################## -# Final processing -############################################## -if [ -e "$pgmout" ] ; then - cat $pgmout -fi - - -msg="ENDED NORMALLY." -postmsg "$jlogfile" "$msg" - - -########################################## -# Remove the Temporary working directory -########################################## -cd $DATAROOT -[[ $KEEPDATA = "NO" ]] && rm -rf $DATA - -date -exit 0 diff --git a/jobs/JGDAS_EFSOI_SFC b/jobs/JGDAS_EFSOI_SFC deleted file mode 100755 index 5924393c33..0000000000 --- a/jobs/JGDAS_EFSOI_SFC +++ /dev/null @@ -1,155 +0,0 @@ -#!/bin/ksh -set -x - -export RUN_ENVIR=${RUN_ENVIR:-"nco"} -export PS4='$SECONDS + ' -date - - -############################# -# Source relevant config files -############################# -export EXPDIR=${EXPDIR:-$HOMEgfs/parm/config} -configs="base esfcfsoi" -config_path=${EXPDIR:-$NWROOT/gfs.${gfs_ver}/parm/config} -for config in $configs; do - . $config_path/config.$config - status=$? - [[ $status -ne 0 ]] && exit $status -done - - -########################################## -# Source machine runtime environment -########################################## -. $HOMEgfs/env/${machine}.env esfcfsoi -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# Obtain unique process id (pid) and make temp directory -############################################## -export pid=${pid:-$$} -export outid=${outid:-"LL$job"} -export DATA=${DATA:-${DATAROOT}/${jobid:?}} -mkdir -p $DATA -cd $DATA - - -############################################## -# Run setpdy and initialize PDY variables -############################################## -export cycle="t${cyc}z" -setpdy.sh -. ./PDY - - -############################################## -# Determine Job Output Name on System -############################################## -export pgmout="OUTPUT.${pid}" -export pgmerr=errfile - - -############################################## -# Set variables used in the script -############################################## -export CDATE=${CDATE:-${PDY}${cyc}} -export CDUMP=${CDUMP:-${RUN:-"gdas"}} -export COMPONENT=${COMPONENT:-atmos} -if [ $RUN_ENVIR = "nco" ]; then - export ROTDIR=${COMROOT:?}/$NET/$envir -fi - - -############################################## -# Begin JOB SPECIFIC work -############################################## - -GDATE=$($NDATE -$assim_freq $CDATE) -gPDY=$(echo $GDATE | cut -c1-8) -gcyc=$(echo $GDATE | cut -c9-10) -GDUMP=${GDUMP:-"gdas"} - -export CASE=$CASE_ENKF - - -EUPD_CYC=$(echo ${EUPD_CYC:-"gdas"} | tr a-z A-Z) -if [ $EUPD_CYC = "GFS" ]; then - CDUMP_ENKF="gfs" -else - CDUMP_ENKF=$CDUMP -fi - -export OPREFIX="${CDUMP}.t${cyc}z." -export APREFIX="${CDUMP}.t${cyc}z." -export APREFIX_ENKF="${CDUMP_ENKF}.t${cyc}z." -export GPREFIX="${CDUMP}.t${gcyc}z." -export GSUFFIX=${GSUFFIX:-$SUFFIX} -export ASUFFIX=${ASUFFIX:-$SUFFIX} - -if [ $RUN_ENVIR = "nco" -o ${ROTDIR_DUMP:-NO} = "YES" ]; then - export COMIN_OBS=${COMIN_OBS:-$ROTDIR/$RUN.$PDY/$cyc/$COMPONENT} - export COMIN_GES_OBS=${COMIN_GES_OBS:-$ROTDIR/$GDUMP.$gPDY/$gcyc/$COMPONENT} -else - export COMIN_OBS="$DMPDIR/$CDATE/$CDUMP" - export COMIN_GES_OBS="$DMPDIR/$GDATE/$GDUMP" -fi - -# COMIN, COMIN_ENS and COMIN_GES_ENS are used in script -# Current GDAS comrot folder -export COMIN="$ROTDIR/$CDUMP.$PDY/$cyc/$COMPONENT" - -# Current ENKF (EFSOI) comrot folder -export COMIN_ENS="$ROTDIR/efsoi$CDUMP_ENKF.$PDY/$cyc/$COMPONENT" - -# Current ENKF (EFSOI) comrot folder -export COMOUT_ENS="$ROTDIR/efsoi$CDUMP.$PDY/$cyc/$COMPONENT" - -# Previous ENKF comrot folder -export COMIN_GES_ENS="$ROTDIR/enkf$CDUMP.$gPDY/$gcyc/$COMPONENT" - -#export COMIN="$ROTDIR/$CDUMP.$PDY/$cyc/$COMPONENT" -#export COMIN_ENS="$ROTDIR/enkf$CDUMP_ENKF.$PDY/$cyc/$COMPONENT" -#export COMOUT_ENS="$ROTDIR/enkf$CDUMP.$PDY/$cyc/$COMPONENT" -#export COMIN_GES_ENS="$ROTDIR/enkf$CDUMP.$gPDY/$gcyc/$COMPONENT" - - -############################################################### -# Run relevant script -env -msg="HAS BEGUN on `hostname`" -postmsg "$jlogfile" "$msg" -$LOGSCRIPT - - -${ENKFRESFCSH:-$SCRgfs/exgdas_enkf_sfc.sh} -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# End JOB SPECIFIC work -############################################## - -############################################## -# Final processing -############################################## -if [ -e "$pgmout" ] ; then - cat $pgmout -fi - - -msg="ENDED NORMALLY." -postmsg "$jlogfile" "$msg" - - -########################################## -# Remove the Temporary working directory -########################################## -cd $DATAROOT -[[ $KEEPDATA = "NO" ]] && rm -rf $DATA - -date -exit 0 diff --git a/jobs/JGDAS_EFSOI_UPDATE b/jobs/JGDAS_EFSOI_UPDATE deleted file mode 100755 index 3d730eab7e..0000000000 --- a/jobs/JGDAS_EFSOI_UPDATE +++ /dev/null @@ -1,137 +0,0 @@ -#!/bin/ksh -set -x - -export RUN_ENVIR=${RUN_ENVIR:-"nco"} -export PS4='$SECONDS + ' -date - - -############################# -# Source relevant config files -############################# -export EXPDIR=${EXPDIR:-$HOMEgfs/parm/config} - -configs="base anal eupdfsoi" -config_path=${EXPDIR:-$NWROOT/gfs.${gfs_ver}/parm/config} -for config in $configs; do - . $config_path/config.$config - status=$? - [[ $status -ne 0 ]] && exit $status -done - - -########################################## -# Source machine runtime environment -########################################## -. $HOMEgfs/env/${machine}.env eupdfsoi - - -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# Obtain unique process id (pid) and make temp directory -############################################## -export pid=${pid:-$$} -export outid=${outid:-"LL$job"} -export DATA=${DATA:-${DATAROOT}/${jobid:?}} -mkdir -p $DATA -cd $DATA - - -############################################## -# Run setpdy and initialize PDY variables -############################################## -export cycle="t${cyc}z" -setpdy.sh -. ./PDY - - -############################################## -# Determine Job Output Name on System -############################################## -export pgmout="OUTPUT.${pid}" -export pgmerr=errfile - - -############################################## -# Set variables used in the exglobal script -############################################## -export CDATE=${CDATE:-${PDY}${cyc}} -export CDUMP=${CDUMP:-${RUN:-"gdas"}} -export COMPONENT=${COMPONENT:-atmos} -if [ $RUN_ENVIR = "nco" ]; then - export ROTDIR=${COMROOT:?}/$NET/$envir -fi - - -############################################## -# Begin JOB SPECIFIC work -############################################## - -GDATE=$($NDATE -$assim_freq $CDATE) -gPDY=$(echo $GDATE | cut -c1-8) -gcyc=$(echo $GDATE | cut -c9-10) - -export APREFIX="${CDUMP}.t${cyc}z." -export GPREFIX="gdas.t${gcyc}z." -export ASUFFIX=${ASUFFIX:-$SUFFIX} -export GSUFFIX=${GSUFFIX:-$SUFFIX} - - - -# COMIN_GES_ENS and COMOUT_ANL_ENS are used in script -export COMIN_GES_ENS="$ROTDIR/enkfgdas.$gPDY/$gcyc/$COMPONENT" -export COMOUT_ANL_ENS="$ROTDIR/enkf$CDUMP.$PDY/$cyc/$COMPONENT" -export COMOUT_ANL_ENSFSOI="$ROTDIR/efsoi$CDUMP.$PDY/$cyc/$COMPONENT" - - - -mkdir -p $COMOUT_ANL_ENSFSOI - -############################################################### -# Run relevant exglobal script -env -msg="HAS BEGUN on `hostname`" -postmsg "$jlogfile" "$msg" -$LOGSCRIPT - -${EFSOIUPDSH:-$SCRgfs/exgdas_efsoi_update.sh} -status=$? -[[ $status -ne 0 ]] && exit $status - - -############################################## -# Send Alerts -############################################## -if [ $SENDDBN = YES ] ; then - $DBNROOT/bin/dbn_alert MODEL ENKF1_MSC_enkfstat $job $COMOUT_ANL_ENSFSOI/${APREFIX}enkfstat -fi - - -############################################## -# End JOB SPECIFIC work -############################################## - -############################################## -# Final processing -############################################## -if [ -e "$pgmout" ] ; then - cat $pgmout -fi - - -msg="ENDED NORMALLY." -postmsg "$jlogfile" "$msg" - - -########################################## -# Remove the Temporary working directory -########################################## -cd $DATAROOT - -[[ $KEEPDATA = "NO" ]] && rm -rf $DATA - -date -exit 0 diff --git a/jobs/JGDAS_ENKF_FCST b/jobs/JGDAS_ENKF_FCST index adc0b492e6..60da5e7063 100755 --- a/jobs/JGDAS_ENKF_FCST +++ b/jobs/JGDAS_ENKF_FCST @@ -110,6 +110,7 @@ if [ $npass -ne $NMEM_EFCSGRP ]; then exit 99 fi + ############################################## # Send Alerts ############################################## diff --git a/scripts/exgdas_efsoi.sh b/scripts/exgdas_efsoi.sh deleted file mode 100755 index 506bba099a..0000000000 --- a/scripts/exgdas_efsoi.sh +++ /dev/null @@ -1,475 +0,0 @@ -#!/bin/ksh -################################################################################ -#### UNIX Script Documentation Block -# . . -# Script name: exgdas_efsoi.sh -# Script description: Make global_enkf update for efsoi -# -# Author: Rahul Mahajan Org: NCEP/EMC Date: 2017-03-02 -# Author: Liaofan Lin/Andrew Eichmann Date: 2020-12-03 -# -# Abstract: This script runs the efsoi executable -# -# $Id$ -# -# Attributes: -# Language: POSIX shell -# Machine: WCOSS-Cray/Theia -# -################################################################################ - -# Set environment. -VERBOSE=${VERBOSE:-"YES"} -if [ $VERBOSE = "YES" ] ; then - echo $(date) EXECUTING $0 $* >&2 - set -x -fi - -# Directories. -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"} -CFP_MP=${CFP_MP:-"NO"} -nm="" -if [ $CFP_MP = "YES" ]; then - nm=0 -fi -APRUNCFP=${APRUNCFP:-""} -APRUN_ENKF=${APRUN_ENKF:-${APRUN:-""}} -NTHREADS_ENKF=${NTHREADS_ENKF:-${NTHREADS:-1}} - -# Executables -#ENKFEXEC=${ENKFEXEC:-$HOMEgfs/exec/global_enkf.x} -EFSOIEXEC=${EFSOIEXEC:-$HOMEgfs/exec/efsoi.x} - -# Cycling and forecast hour specific parameters -CDATE=${CDATE:-"2001010100"} - -# Filenames. -GPREFIX=${GPREFIX:-""} -GSUFFIX=${GSUFFIX:-$SUFFIX} -APREFIX=${APREFIX:-""} -ASUFFIX=${ASUFFIX:-$SUFFIX} - -SMOOTH_ENKF=${SMOOTH_ENKF:-"YES"} - -GBIASe=${GBIASe:-${APREFIX}abias_int.ensmean} -CNVSTAT=${CNVSTAT:-${APREFIX}cnvstat} # AFE not needed? -OZNSTAT=${OZNSTAT:-${APREFIX}oznstat} # AFE not needed? -RADSTAT=${RADSTAT:-${APREFIX}radstat} # AFE not needed? -#ENKFSTAT=${ENKFSTAT:-${APREFIX}enkfstat} -EFSOISTAT=${EFSOISTAT:-${APREFIX}efsoistat} - -#AFE for EFSOI -#VERFANL=${VERFANL:-${APREFIX}atmanl.nemsio} -VERFANL=${VERFANL:-${APREFIX}atmanl.ensres.nc} -FCSTLONG=${GPREFIX}atmf030.ensmean.nc -FCSTSHORT=${APREFIX}atmf024.ensmean.nc -FCST6HRS=${PPREFIX}atmf006.ensmean.nc -OSENSEIN=osense_${CDATE}_init.dat -OSENSEOUT=osense_${CDATE}_final.dat - -# this needs to be set manually because params in enkf will default to fhr03 -fgfileprefixes=sfg_${CDATE}_fhr06_ - - - - -#analysise Namelst parameters -USE_CORRELATED_OBERRS=${USE_CORRELATED_OBERRS:-"NO"} -NMEM_ENKF=${NMEM_ENKF:-80} -NAM_ENKF=${NAM_ENKF:-""} -SATOBS_ENKF=${SATOBS_ENKF:-""} -OZOBS_ENKF=${OZOBS_ENKF:-""} -use_correlated_oberrs=${use_correlated_oberrs:-".false."} -if [ $USE_CORRELATED_OBERRS == "YES" ]; then - use_correlated_oberrs=".true." -fi -imp_physics=${imp_physics:-"99"} -#lupp=${lupp:-".true."} -lupp=${lupp:-".false."} # AFE to match old EFSOI -corrlength=${corrlength:-1250} -lnsigcutoff=${lnsigcutoff:-2.5} -analpertwt=${analpertwt:-0.85} -#readin_localization_enkf=${readin_localization_enkf:-".true."} -readin_localization_enkf=${readin_localization_enkf:-".false."} # AFE -reducedgrid=${reducedgrid:-".true."} -letkf_flag=${letkf_flag:-".false."} -getkf=${getkf:-".false."} -denkf=${denkf:-".false."} -nobsl_max=${nobsl_max:-10000} -lobsdiag_forenkf=${lobsdiag_forenkf:-".false."} -write_spread_diag=${write_spread_diag:-".false."} -cnvw_option=${cnvw_option:-".false."} -netcdf_diag=${netcdf_diag:-".true."} -modelspace_vloc=${modelspace_vloc:-".false."} # if true, 'vlocal_eig.dat' is needed -IAUFHRS_ENKF=${IAUFHRS_ENKF:-6} -DO_CALC_INCREMENT=${DO_CALC_INCREMENT:-"NO"} -INCREMENTS_TO_ZERO=${INCREMENTS_TO_ZERO:-"'NONE'"} - -################################################################################ -#ATMGES_ENSMEAN=$COMIN_GES_ENS/${GPREFIX}atmf006.ensmean${GSUFFIX} -ATMGES_ENSMEAN=$COMIN_ANL/$VERFANL -if [ $SUFFIX = ".nc" ]; then - LONB_ENKF=${LONB_ENKF:-$($NCLEN $ATMGES_ENSMEAN grid_xt)} # get LONB_ENKF - LATB_ENKF=${LATB_ENKF:-$($NCLEN $ATMGES_ENSMEAN grid_yt)} # get LATB_ENFK - LEVS_ENKF=${LEVS_ENKF:-$($NCLEN $ATMGES_ENSMEAN pfull)} # get LEVS_ENFK - use_gfs_ncio=".true." - use_gfs_nemsio=".false." - paranc=${paranc:-".true."} - if [ $DO_CALC_INCREMENT = "YES" ]; then - write_fv3_incr=".false." - else - write_fv3_incr=".true." - WRITE_INCR_ZERO="incvars_to_zero= $INCREMENTS_TO_ZERO," - fi -else - LEVS_ENKF=${LEVS_ENKF:-$($NEMSIOGET $ATMGES_ENSMEAN dimz | awk '{print $2}')} - LATB_ENKF=${LATB_ENKF:-$($NEMSIOGET $ATMGES_ENSMEAN dimy | awk '{print $2}')} - LONB_ENKF=${LONB_ENKF:-$($NEMSIOGET $ATMGES_ENSMEAN dimx | awk '{print $2}')} - use_gfs_ncio=".false." - use_gfs_nemsio=".true." - paranc=${paranc:-".false."} -fi -LATA_ENKF=${LATA_ENKF:-$LATB_ENKF} -LONA_ENKF=${LONA_ENKF:-$LONB_ENKF} - -SATANGL=${SATANGL:-${FIXgsi}/global_satangbias.txt} -SATINFO=${SATINFO:-${FIXgsi}/global_satinfo.txt} -CONVINFO=${CONVINFO:-${FIXgsi}/global_convinfo.txt} -OZINFO=${OZINFO:-${FIXgsi}/global_ozinfo.txt} -SCANINFO=${SCANINFO:-${FIXgsi}/global_scaninfo.txt} -HYBENSINFO=${HYBENSINFO:-${FIXgsi}/global_hybens_info.l${LEVS_ENKF}.txt} -ANAVINFO=${ANAVINFO:-${FIXgsi}/global_anavinfo.l${LEVS_ENKF}.txt} -VLOCALEIG=${VLOCALEIG:-${FIXgsi}/vlocal_eig_l${LEVS_ENKF}.dat} - -ENKF_SUFFIX="s" -[[ $SMOOTH_ENKF = "NO" ]] && ENKF_SUFFIX="" - -################################################################################ -# Preprocessing -mkdata=NO -if [ ! -d $DATA ]; then - mkdata=YES - mkdir -p $DATA -fi -cd $DATA || exit 99 - -################################################################################ -# Fixed files -$NLN $SATANGL satbias_angle -$NLN $SATINFO satinfo -$NLN $SCANINFO scaninfo -$NLN $CONVINFO convinfo -$NLN $OZINFO ozinfo -$NLN $HYBENSINFO hybens_info -$NLN $ANAVINFO anavinfo -$NLN $VLOCALEIG vlocal_eig.dat - -# Bias correction coefficients based on the ensemble mean -$NLN $COMOUT_ANL_ENSFSOI/$GBIASe satbias_in - -################################################################################ - -#if [ $USE_CFP = "YES" ]; then -# [[ -f $DATA/untar.sh ]] && rm $DATA/untar.sh -# [[ -f $DATA/mp_untar.sh ]] && rm $DATA/mp_untar.sh -# set +x -# cat > $DATA/untar.sh << EOFuntar -##!/bin/sh -#memchar=\$1 -#flist="$CNVSTAT $OZNSTAT $RADSTAT" -#for ftype in \$flist; do -# if [ \$memchar = "ensmean" ]; then -# fname=$COMOUT_ANL_ENS/\${ftype}.ensmean -# else -# fname=$COMOUT_ANL_ENS/\$memchar/\$ftype -# fi -# tar -xvf \$fname -#done -#EOFuntar -# set -x -# chmod 755 $DATA/untar.sh -#fi - -################################################################################ -# Ensemble guess, observational data and analyses/increments - -#flist="$CNVSTAT $OZNSTAT $RADSTAT" -#if [ $USE_CFP = "YES" ]; then -# echo "$nm $DATA/untar.sh ensmean" | tee -a $DATA/mp_untar.sh -# if [ ${CFP_MP:-"NO"} = "YES" ]; then -# nm=$((nm+1)) -# fi -#else -# for ftype in $flist; do -# fname=$COMOUT_ANL_ENS/${ftype}.ensmean -# tar -xvf $fname -# done -#fi -nfhrs=`echo $IAUFHRS_ENKF | sed 's/,/ /g'` -for imem in $(seq 1 $NMEM_ENKF); do - memchar="mem"$(printf %03i $imem) - mkdir ${memchar} - $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atmf024.nc ${memchar} -# if [ $lobsdiag_forenkf = ".false." ]; then -# if [ $USE_CFP = "YES" ]; then -# echo "$nm $DATA/untar.sh $memchar" | tee -a $DATA/mp_untar.sh -# if [ ${CFP_MP:-"NO"} = "YES" ]; then -# nm=$((nm+1)) -# fi -# else -# for ftype in $flist; do -# fname=$COMOUT_ANL_ENS/$memchar/$ftype -# tar -xvf $fname -# done -# fi -# fi -##### mkdir -p $COMOUT_ANL_ENSFSOI/$memchar -# for FHR in $nfhrs; do -# $NLN $COMIN_GES_ENS/$memchar/${GPREFIX}atmf00${FHR}${ENKF_SUFFIX}${GSUFFIX} sfg_${CDATE}_fhr0${FHR}_${memchar} -# if [ $cnvw_option = ".true." ]; then -# $NLN $COMIN_GES_ENS/$memchar/${GPREFIX}sfcf00${FHR}${GSUFFIX} sfgsfc_${CDATE}_fhr0${FHR}_${memchar} -# fi -# if [ $FHR -eq 6 ]; then -# if [ $DO_CALC_INCREMENT = "YES" ]; then -# #$NLN $COMOUT_ANL_ENS/$memchar/${APREFIX}atmanl${ASUFFIX} sanl_${CDATE}_fhr0${FHR}_${memchar} -# $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atmanl${ASUFFIX} sanl_${CDATE}_fhr0${FHR}_ni${memchar} -# else -# #$NLN $COMOUT_ANL_ENS/$memchar/${APREFIX}atminc${ASUFFIX} incr_${CDATE}_fhr0${FHR}_${memchar} -# $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atminc${ASUFFIX} incr_${CDATE}_fhr0${FHR}_ni${memchar} -# fi -# else -# if [ $DO_CALC_INCREMENT = "YES" ]; then -# #$NLN $COMOUT_ANL_ENS/$memchar/${APREFIX}atma00${FHR}${ASUFFIX} sanl_${CDATE}_fhr0${FHR}_${memchar} -# $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atma00${FHR}${ASUFFIX} sanl_${CDATE}_fhr0${FHR}_ni${memchar} -# else -# #$NLN $COMOUT_ANL_ENS/$memchar/${APREFIX}atmi00${FHR}${ASUFFIX} incr_${CDATE}_fhr0${FHR}_${memchar} -# $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atmi00${FHR}${ASUFFIX} incr_${CDATE}_fhr0${FHR}_ni${memchar} -# fi -# fi -# done -done - -# Ensemble mean guess -#for FHR in $nfhrs; do -# $NLN $COMIN_GES_ENS/${GPREFIX}atmf00${FHR}.ensmean${GSUFFIX} sfg_${CDATE}_fhr0${FHR}_ensmean -# if [ $cnvw_option = ".true." ]; then -# $NLN $COMIN_GES_ENS/${GPREFIX}sfcf00${FHR}.ensmean${GSUFFIX} sfgsfc_${CDATE}_fhr0${FHR}_ensmean -# fi -#done - -$NLN $COMIN_GES_ENS/${GPREFIX}atmf006.ensmean${GSUFFIX} sfg_${CDATE}_fhr06_ensmean -$NLN $COMIN_GES_ENS/${GPREFIX}atmf006.ensmean${GSUFFIX} sfg_${CDATE}_fhr03_ensmean - -# verifying analysis -$NLN ${COMIN_ANL}/${VERFANL} . -# the above way is the proper one, but this bit of subterfuge -# is necessary because ensmean is hard coded into GSI - needs -# to changed -$NLN ${COMIN_ANL}/${VERFANL} ${APREFIX}atmanl.ensmean.nc - -# forecasts -$NLN $COMIN_GES_ENS/$FCSTLONG . -$NLN $COMOUT_ANL_ENSFSOI/$FCSTSHORT . - -# inital osense file -# efsoi.x will read then clobber this -$NCP $COMOUT_ANL_ENSFSOI/$OSENSEIN osense_${CDATE}.dat - - - - -if [ $USE_CFP = "YES" ]; then - chmod 755 $DATA/mp_untar.sh - ncmd=$(cat $DATA/mp_untar.sh | wc -l) - if [ $ncmd -gt 0 ]; 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 - fi -fi - -################################################################################ -# Create global_enkf namelist -# AFE changed from original: -# gdatehr, datehr, andataname added -# analpertwt and lnsigcutoff changed upstream -# numiter from 0 to 1 -# fso_cycling=.true., -# efsoi_flag=.true., -# wmoist=1.0,adrate=0.75 - - - -cat > enkf.nml << EOFnml -&nam_enkf - datestring="$CDATE",datapath="$DATA/", - gdatehr=$gcyc, - datehr=$cyc, - fgfileprefixes=$fgfileprefixes - andataname="$VERFANL", - analpertwtnh=${analpertwt},analpertwtsh=${analpertwt},analpertwttr=${analpertwt}, - covinflatemax=1.e2,covinflatemin=1,pseudo_rh=.true.,iassim_order=0, - corrlengthnh=${corrlength},corrlengthsh=${corrlength},corrlengthtr=${corrlength}, - lnsigcutoffnh=${lnsigcutoff},lnsigcutoffsh=${lnsigcutoff},lnsigcutofftr=${lnsigcutoff}, - lnsigcutoffpsnh=${lnsigcutoff},lnsigcutoffpssh=${lnsigcutoff},lnsigcutoffpstr=${lnsigcutoff}, - lnsigcutoffsatnh=${lnsigcutoff},lnsigcutoffsatsh=${lnsigcutoff},lnsigcutoffsattr=${lnsigcutoff}, - obtimelnh=1.e30,obtimelsh=1.e30,obtimeltr=1.e30, - saterrfact=1.0,numiter=1, - sprd_tol=1.e30,paoverpb_thresh=0.98, - nlons=$LONA_ENKF,nlats=$LATA_ENKF,nlevs=$LEVS_ENKF,nanals=$NMEM_ENKF, - deterministic=.true.,sortinc=.true.,lupd_satbiasc=.false., - reducedgrid=${reducedgrid},readin_localization=${readin_localization_enkf}., - use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},imp_physics=$imp_physics,lupp=$lupp, - univaroz=.false.,adp_anglebc=.true.,angord=4,use_edges=.false.,emiss_bc=.true., - letkf_flag=${letkf_flag},nobsl_max=${nobsl_max},denkf=${denkf},getkf=${getkf}., - nhr_anal=${IAUFHRS_ENKF},nhr_state=${IAUFHRS_ENKF},use_qsatensmean=.true., - lobsdiag_forenkf=$lobsdiag_forenkf, - write_spread_diag=$write_spread_diag, - modelspace_vloc=$modelspace_vloc, - use_correlated_oberrs=${use_correlated_oberrs}, - netcdf_diag=$netcdf_diag,cnvw_option=$cnvw_option, - paranc=$paranc,write_fv3_incr=$write_fv3_incr, - efsoi_cycling=.true., - efsoi_flag=.true., - wmoist=1.0,adrate=0.75 - $WRITE_INCR_ZERO - $NAM_ENKF -/ -&satobs_enkf - sattypes_rad(1) = 'amsua_n15', dsis(1) = 'amsua_n15', - sattypes_rad(2) = 'amsua_n18', dsis(2) = 'amsua_n18', - sattypes_rad(3) = 'amsua_n19', dsis(3) = 'amsua_n19', - sattypes_rad(4) = 'amsub_n16', dsis(4) = 'amsub_n16', - sattypes_rad(5) = 'amsub_n17', dsis(5) = 'amsub_n17', - sattypes_rad(6) = 'amsua_aqua', dsis(6) = 'amsua_aqua', - sattypes_rad(7) = 'amsua_metop-a', dsis(7) = 'amsua_metop-a', - sattypes_rad(8) = 'airs_aqua', dsis(8) = 'airs_aqua', - sattypes_rad(9) = 'hirs3_n17', dsis(9) = 'hirs3_n17', - sattypes_rad(10)= 'hirs4_n19', dsis(10)= 'hirs4_n19', - sattypes_rad(11)= 'hirs4_metop-a', dsis(11)= 'hirs4_metop-a', - sattypes_rad(12)= 'mhs_n18', dsis(12)= 'mhs_n18', - sattypes_rad(13)= 'mhs_n19', dsis(13)= 'mhs_n19', - sattypes_rad(14)= 'mhs_metop-a', dsis(14)= 'mhs_metop-a', - sattypes_rad(15)= 'goes_img_g11', dsis(15)= 'imgr_g11', - sattypes_rad(16)= 'goes_img_g12', dsis(16)= 'imgr_g12', - sattypes_rad(17)= 'goes_img_g13', dsis(17)= 'imgr_g13', - sattypes_rad(18)= 'goes_img_g14', dsis(18)= 'imgr_g14', - sattypes_rad(19)= 'goes_img_g15', dsis(19)= 'imgr_g15', - sattypes_rad(20)= 'avhrr_n18', dsis(20)= 'avhrr3_n18', - sattypes_rad(21)= 'avhrr_metop-a', dsis(21)= 'avhrr3_metop-a', - sattypes_rad(22)= 'avhrr_n19', dsis(22)= 'avhrr3_n19', - sattypes_rad(23)= 'amsre_aqua', dsis(23)= 'amsre_aqua', - sattypes_rad(24)= 'ssmis_f16', dsis(24)= 'ssmis_f16', - sattypes_rad(25)= 'ssmis_f17', dsis(25)= 'ssmis_f17', - sattypes_rad(26)= 'ssmis_f18', dsis(26)= 'ssmis_f18', - sattypes_rad(27)= 'ssmis_f19', dsis(27)= 'ssmis_f19', - sattypes_rad(28)= 'ssmis_f20', dsis(28)= 'ssmis_f20', - sattypes_rad(29)= 'sndrd1_g11', dsis(29)= 'sndrD1_g11', - sattypes_rad(30)= 'sndrd2_g11', dsis(30)= 'sndrD2_g11', - sattypes_rad(31)= 'sndrd3_g11', dsis(31)= 'sndrD3_g11', - sattypes_rad(32)= 'sndrd4_g11', dsis(32)= 'sndrD4_g11', - sattypes_rad(33)= 'sndrd1_g12', dsis(33)= 'sndrD1_g12', - sattypes_rad(34)= 'sndrd2_g12', dsis(34)= 'sndrD2_g12', - sattypes_rad(35)= 'sndrd3_g12', dsis(35)= 'sndrD3_g12', - sattypes_rad(36)= 'sndrd4_g12', dsis(36)= 'sndrD4_g12', - sattypes_rad(37)= 'sndrd1_g13', dsis(37)= 'sndrD1_g13', - sattypes_rad(38)= 'sndrd2_g13', dsis(38)= 'sndrD2_g13', - sattypes_rad(39)= 'sndrd3_g13', dsis(39)= 'sndrD3_g13', - sattypes_rad(40)= 'sndrd4_g13', dsis(40)= 'sndrD4_g13', - sattypes_rad(41)= 'sndrd1_g14', dsis(41)= 'sndrD1_g14', - sattypes_rad(42)= 'sndrd2_g14', dsis(42)= 'sndrD2_g14', - sattypes_rad(43)= 'sndrd3_g14', dsis(43)= 'sndrD3_g14', - sattypes_rad(44)= 'sndrd4_g14', dsis(44)= 'sndrD4_g14', - sattypes_rad(45)= 'sndrd1_g15', dsis(45)= 'sndrD1_g15', - sattypes_rad(46)= 'sndrd2_g15', dsis(46)= 'sndrD2_g15', - sattypes_rad(47)= 'sndrd3_g15', dsis(47)= 'sndrD3_g15', - sattypes_rad(48)= 'sndrd4_g15', dsis(48)= 'sndrD4_g15', - sattypes_rad(49)= 'iasi_metop-a', dsis(49)= 'iasi_metop-a', - sattypes_rad(50)= 'seviri_m08', dsis(50)= 'seviri_m08', - sattypes_rad(51)= 'seviri_m09', dsis(51)= 'seviri_m09', - sattypes_rad(52)= 'seviri_m10', dsis(52)= 'seviri_m10', - sattypes_rad(53)= 'seviri_m11', dsis(53)= 'seviri_m11', - sattypes_rad(54)= 'amsua_metop-b', dsis(54)= 'amsua_metop-b', - sattypes_rad(55)= 'hirs4_metop-b', dsis(55)= 'hirs4_metop-b', - sattypes_rad(56)= 'mhs_metop-b', dsis(56)= 'mhs_metop-b', - sattypes_rad(57)= 'iasi_metop-b', dsis(57)= 'iasi_metop-b', - sattypes_rad(58)= 'avhrr_metop-b', dsis(58)= 'avhrr3_metop-b', - sattypes_rad(59)= 'atms_npp', dsis(59)= 'atms_npp', - sattypes_rad(60)= 'atms_n20', dsis(60)= 'atms_n20', - sattypes_rad(61)= 'cris_npp', dsis(61)= 'cris_npp', - sattypes_rad(62)= 'cris-fsr_npp', dsis(62)= 'cris-fsr_npp', - sattypes_rad(63)= 'cris-fsr_n20', dsis(63)= 'cris-fsr_n20', - sattypes_rad(64)= 'gmi_gpm', dsis(64)= 'gmi_gpm', - sattypes_rad(65)= 'saphir_meghat', dsis(65)= 'saphir_meghat', - sattypes_rad(66)= 'amsua_metop-c', dsis(66)= 'amsua_metop-c', - sattypes_rad(67)= 'mhs_metop-c', dsis(67)= 'mhs_metop-c', - sattypes_rad(68)= 'ahi_himawari8', dsis(68)= 'ahi_himawari8', - sattypes_rad(69)= 'abi_g16', dsis(69)= 'abi_g16', - sattypes_rad(70)= 'abi_g17', dsis(70)= 'abi_g17', - $SATOBS_ENKF -/ -&ozobs_enkf - sattypes_oz(1) = 'sbuv2_n16', - sattypes_oz(2) = 'sbuv2_n17', - sattypes_oz(3) = 'sbuv2_n18', - sattypes_oz(4) = 'sbuv2_n19', - sattypes_oz(5) = 'omi_aura', - sattypes_oz(6) = 'gome_metop-a', - sattypes_oz(7) = 'gome_metop-b', - sattypes_oz(8) = 'mls30_aura', - sattypes_oz(9) = 'ompsnp_npp', - sattypes_oz(10) = 'ompstc8_npp', - $OZOBS_ENKF -/ -EOFnml - -################################################################################ -# Run enkf update - -export OMP_NUM_THREADS=$NTHREADS_ENKF -#export pgm=$ENKFEXEC AFE -export pgm=$EFSOIEXEC -. prep_step - -#$NCP $ENKFEXEC $DATA AFE -#$APRUN_ENKF ${DATA}/$(basename $ENKFEXEC) 1>stdout 2>stderr AFE -$NCP $EFSOIEXEC $DATA -$APRUN_ENKF ${DATA}/$(basename $EFSOIEXEC) 1>stdout 2>stderr -rc=$? - -export ERR=$rc -export err=$ERR -$ERRSCRIPT || exit 2 - -## save for EFSOI task (still needed?) -#$NCP $COMOUT_ANL_ENS/$GBIASe $COMOUT_ANL_ENSFSOI - -# Cat runtime output files. -cat stdout stderr > $COMOUT_ANL_ENSFSOI/$EFSOISTAT - -$NCP osense_${CDATE}.dat $COMOUT_ANL_ENSFSOI/$OSENSEOUT -$NCP osense_${CDATE}.dat $OSENSE_SAVE_DIR/$OSENSEOUT - -################################################################################ -# Postprocessing -######## AFE remove after testing -cp -r $DATA /scratch1/NCEPDEV/stmp4/Andrew.Eichmann/efsoi -######## AFE -cd $pwd -[[ $mkdata = "YES" ]] && rm -rf $DATA -set +x -if [ $VERBOSE = "YES" ]; then - echo $(date) EXITING $0 with return code $err >&2 -fi -exit $err diff --git a/scripts/exgdas_efsoi_fcst.sh b/scripts/exgdas_efsoi_fcst.sh deleted file mode 100755 index 13f75cd625..0000000000 --- a/scripts/exgdas_efsoi_fcst.sh +++ /dev/null @@ -1,243 +0,0 @@ -#!/bin/ksh -################################################################################ -#### UNIX Script Documentation Block -# . . -# Script name: exgdas_efsoi_fcst.sh -# Script description: Run ensemble forecasts -# -# Author: Rahul Mahajan Org: NCEP/EMC Date: 2017-03-02 -# Author: Andrew Eichmann Org: NCEP/EMC Date: 2021-02-25 -# -# Abstract: This script runs ensemble forecasts serially one-after-another -# for EFSOI. Based on exgdas_enkf_fcst.sh -# -# $Id$ -# -# Attributes: -# Language: POSIX shell -# Machine: WCOSS-Cray/Theia -# -#### -################################################################################ - -# Set environment. -export VERBOSE=${VERBOSE:-"YES"} -if [ $VERBOSE = "YES" ] ; then - echo $(date) EXECUTING $0 $* >&2 - set -x -fi - -# Directories. -pwd=$(pwd) -export FIX_DIR=${FIX_DIR:-$HOMEgfs/fix} -export FIX_AM=${FIX_AM:-$FIX_DIR/fix_am} - -# Utilities -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} - -# Enemble group, begin and end -ENSGRP=${ENSGRP:-1} -ENSBEG=${ENSBEG:-1} -ENSEND=${ENSEND:-1} - -# Model builds -export FCSTEXECDIR=${FCSTEXECDIR:-$HOMEgfs/sorc/fv3gfs.fd/BUILD/bin} -export FCSTEXEC=${FCSTEXEC:-fv3gfs.x} - -# Get DA specific diag table. -export PARM_FV3DIAG=${PARM_FV3DIAG:-$HOMEgfs/parm/parm_fv3diag} -export DIAG_TABLE=${DIAG_TABLE_ENKF:-${DIAG_TABLE:-$PARM_FV3DIAG/diag_table_da}} - -# Cycling and forecast hour specific parameters -export CDATE=${CDATE:-"2001010100"} -export CDUMP=${CDUMP:-"gdas"} - -# Re-run failed members, or entire group -RERUN_EFCSGRP=${RERUN_EFCSGRP:-"YES"} - -# Recenter flag and increment file prefix -RECENTER_ENKF=${RECENTER_ENKF:-"YES"} -export PREFIX_ATMINC=${PREFIX_ATMINC:-""} - -# Ops related stuff -SENDECF=${SENDECF:-"NO"} -SENDDBN=${SENDDBN:-"NO"} -GSUFFIX=${GSUFFIX:-$SUFFIX} - -################################################################################ -# Preprocessing -mkdata=NO -if [ ! -d $DATA ]; then - mkdata=YES - mkdir -p $DATA -fi -cd $DATA || exit 99 -DATATOP=$DATA - -################################################################################ -# Set output data -cymd=$(echo $CDATE | cut -c1-8) -chh=$(echo $CDATE | cut -c9-10) -EFCSGRP=$COMOUT/efcs.grp${ENSGRP} -if [ -f $EFCSGRP ]; then - if [ $RERUN_EFCSGRP = "YES" ]; then - rm -f $EFCSGRP - else - echo "RERUN_EFCSGRP = $RERUN_EFCSGRP, will re-run FAILED members only!" - $NMV $EFCSGRP ${EFCSGRP}.fail - fi -fi - -################################################################################ -# Set namelist/model config options common to all members once - -# There are many many model namelist options -# Some are resolution (CASE) dependent, some depend on the model configuration -# and will need to be added here before $FORECASTSH is called -# For now assume that -# 1. the ensemble and the deterministic are same resolution -# 2. the ensemble runs with the same configuration as the deterministic - -# Model config option for Ensemble -export TYPE=${TYPE_ENKF:-${TYPE:-nh}} # choices: nh, hydro -export MONO=${MONO_ENKF:-${MONO:-non-mono}} # choices: mono, non-mono - -# fv_core_nml -export CASE=${CASE_ENKF:-${CASE:-C768}} -export layout_x=${layout_x_ENKF:-${layout_x:-8}} -export layout_y=${layout_y_ENKF:-${layout_y:-16}} -export LEVS=${LEVS_ENKF:-${LEVS:-64}} - -# nggps_diag_nml -export FHOUT=${FHOUT_ENKF:-3} - -# model_configure -export DELTIM=${DELTIM_ENKF:-${DELTIM:-225}} -export FHMAX=${FHMAX_ENKF:-9} -export restart_interval=${restart_interval_ENKF:-${restart_interval:-6}} - -# gfs_physics_nml -export FHSWR=${FHSWR_ENKF:-${FHSWR:-3600.}} -export FHLWR=${FHLWR_ENKF:-${FHLWR:-3600.}} -export IEMS=${IEMS_ENKF:-${IEMS:-1}} -export ISOL=${ISOL_ENKF:-${ISOL:-2}} -export IAER=${IAER_ENKF:-${IAER:-111}} -export ICO2=${ICO2_ENKF:-${ICO2:-2}} -export cdmbgwd=${cdmbgwd_ENKF:-${cdmbgwd:-"3.5,0.25"}} -export dspheat=${dspheat_ENKF:-${dspheat:-".true."}} -export shal_cnv=${shal_cnv_ENKF:-${shal_cnv:-".true."}} -export FHZER=${FHZER_ENKF:-${FHZER:-6}} -export FHCYC=${FHCYC_ENKF:-${FHCYC:-6}} - -# Set PREFIX_ATMINC to r when recentering on -if [ $RECENTER_ENKF = "YES" ]; then - export PREFIX_ATMINC="r" -fi - -# APRUN for different executables -export APRUN_FV3=${APRUN_FV3:-${APRUN:-""}} -export NTHREADS_FV3=${NTHREADS_FV3:-${NTHREADS:-1}} - -################################################################################ -# Run forecast for ensemble member -rc=0 -for imem in $(seq $ENSBEG $ENSEND); do - - cd $DATATOP - - cmem=$(printf %03i $imem) - memchar="mem$cmem" - - echo "Processing MEMBER: $cmem" - - ra=0 - - skip_mem="NO" - if [ -f ${EFCSGRP}.fail ]; then - memstat=$(cat ${EFCSGRP}.fail | grep "MEMBER $cmem" | grep "PASS" | wc -l) - [[ $memstat -eq 1 ]] && skip_mem="YES" - fi - - if [ $skip_mem = "NO" ]; then - - ra=0 - - export MEMBER=$imem - export DATA=$DATATOP/$memchar - if [ -d $DATA ]; then rm -rf $DATA; fi - $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 - fi - - ((rc+=ra)) - - fi - - if [ $SENDDBN = YES ]; then - fhr=$FHOUT - while [ $fhr -le $FHMAX ]; do - FH3=$(printf %03i $fhr) - if [ $(expr $fhr % 3) -eq 0 ]; then - $DBNROOT/bin/dbn_alert MODEL GFS_ENKF $job $COMOUT/$memchar/${CDUMP}.t${cyc}z.sfcf${FH3}${GSUFFIX} - fi - fhr=$((fhr+FHOUT)) - done - fi - - cd $DATATOP - - if [ -s $EFCSGRP ]; then - $NCP $EFCSGRP log_old - fi - [[ -f log ]] && rm log - [[ -f log_new ]] && rm log_new - if [ $ra -ne 0 ]; then - echo "MEMBER $cmem : FAIL" > log - else - echo "MEMBER $cmem : PASS" > log - fi - if [ -s log_old ] ; then - cat log_old log > log_new - else - cat log > log_new - fi - $NCP log_new $EFCSGRP - -done - -################################################################################ -# Echo status of ensemble group -cd $DATATOP -echo "Status of ensemble members in group $ENSGRP:" -cat $EFCSGRP -[[ -f ${EFCSGRP}.fail ]] && rm ${EFCSGRP}.fail - -################################################################################ -# If any members failed, error out -export ERR=$rc -export err=$ERR -$ERRSCRIPT || exit 2 - -################################################################################ -# Postprocessing -cd $pwd -[[ $mkdata = "YES" ]] && rm -rf $DATATOP -set +x -if [ $VERBOSE = "YES" ] ; then - echo $(date) EXITING $0 with return code $err >&2 -fi -exit $err diff --git a/scripts/exgdas_efsoi_update.sh b/scripts/exgdas_efsoi_update.sh deleted file mode 100755 index 3ad322d3c8..0000000000 --- a/scripts/exgdas_efsoi_update.sh +++ /dev/null @@ -1,413 +0,0 @@ -#!/bin/ksh -################################################################################ -#### UNIX Script Documentation Block -# . . -# Script name: exgdas_efsoi_update.sh -# Script description: Make global_enkf update for efsoi -# -# Author: Rahul Mahajan Org: NCEP/EMC Date: 2017-03-02 -# Author: Liaofan Lin/Andrew Eichmann Date: 2020-12-03 -# -# Abstract: This script runs the global_enkf update for efsoi -# -# $Id$ -# -# Attributes: -# Language: POSIX shell -# Machine: WCOSS-Cray/Theia -# -################################################################################ - -# Set environment. -VERBOSE=${VERBOSE:-"YES"} -if [ $VERBOSE = "YES" ] ; then - echo $(date) EXECUTING $0 $* >&2 - set -x -fi - -# Directories. -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"} -CFP_MP=${CFP_MP:-"NO"} -nm="" -if [ $CFP_MP = "YES" ]; then - nm=0 -fi -APRUNCFP=${APRUNCFP:-""} -APRUN_ENKF=${APRUN_ENKF:-${APRUN:-""}} -NTHREADS_ENKF=${NTHREADS_ENKF:-${NTHREADS:-1}} - -# Executables -ENKFEXEC=${ENKFEXEC:-$HOMEgfs/exec/global_enkf.x} - -# Cycling and forecast hour specific parameters -CDATE=${CDATE:-"2001010100"} - -# Filenames. -GPREFIX=${GPREFIX:-""} -GSUFFIX=${GSUFFIX:-$SUFFIX} -APREFIX=${APREFIX:-""} -ASUFFIX=${ASUFFIX:-$SUFFIX} - -SMOOTH_ENKF=${SMOOTH_ENKF:-"YES"} - -GBIASe=${GBIASe:-${APREFIX}abias_int.ensmean} -CNVSTAT=${CNVSTAT:-${APREFIX}cnvstat} -OZNSTAT=${OZNSTAT:-${APREFIX}oznstat} -RADSTAT=${RADSTAT:-${APREFIX}radstat} -ENKFSTAT=${ENKFSTAT:-${APREFIX}enkfstat} - -# Namelist parameters -USE_CORRELATED_OBERRS=${USE_CORRELATED_OBERRS:-"NO"} -NMEM_ENKF=${NMEM_ENKF:-80} -NAM_ENKF=${NAM_ENKF:-""} -SATOBS_ENKF=${SATOBS_ENKF:-""} -OZOBS_ENKF=${OZOBS_ENKF:-""} -use_correlated_oberrs=${use_correlated_oberrs:-".false."} -if [ $USE_CORRELATED_OBERRS == "YES" ]; then - use_correlated_oberrs=".true." -fi -imp_physics=${imp_physics:-"99"} -lupp=${lupp:-".true."} -corrlength=${corrlength:-1250} -lnsigcutoff=${lnsigcutoff:-2.5} -analpertwt=${analpertwt:-0.85} -readin_localization_enkf=${readin_localization_enkf:-".true."} -reducedgrid=${reducedgrid:-".true."} -letkf_flag=${letkf_flag:-".false."} -getkf=${getkf:-".false."} -denkf=${denkf:-".false."} -nobsl_max=${nobsl_max:-10000} -lobsdiag_forenkf=${lobsdiag_forenkf:-".false."} -write_spread_diag=${write_spread_diag:-".false."} -cnvw_option=${cnvw_option:-".false."} -netcdf_diag=${netcdf_diag:-".true."} -modelspace_vloc=${modelspace_vloc:-".false."} # if true, 'vlocal_eig.dat' is needed -IAUFHRS_ENKF=${IAUFHRS_ENKF:-6} -DO_CALC_INCREMENT=${DO_CALC_INCREMENT:-"NO"} -INCREMENTS_TO_ZERO=${INCREMENTS_TO_ZERO:-"'NONE'"} - -################################################################################ -ATMGES_ENSMEAN=$COMIN_GES_ENS/${GPREFIX}atmf006.ensmean${GSUFFIX} -if [ $SUFFIX = ".nc" ]; then - LONB_ENKF=${LONB_ENKF:-$($NCLEN $ATMGES_ENSMEAN grid_xt)} # get LONB_ENKF - LATB_ENKF=${LATB_ENKF:-$($NCLEN $ATMGES_ENSMEAN grid_yt)} # get LATB_ENFK - LEVS_ENKF=${LEVS_ENKF:-$($NCLEN $ATMGES_ENSMEAN pfull)} # get LEVS_ENFK - use_gfs_ncio=".true." - use_gfs_nemsio=".false." - paranc=${paranc:-".true."} - if [ $DO_CALC_INCREMENT = "YES" ]; then - write_fv3_incr=".false." - else - write_fv3_incr=".true." - WRITE_INCR_ZERO="incvars_to_zero= $INCREMENTS_TO_ZERO," - fi -else - LEVS_ENKF=${LEVS_ENKF:-$($NEMSIOGET $ATMGES_ENSMEAN dimz | awk '{print $2}')} - LATB_ENKF=${LATB_ENKF:-$($NEMSIOGET $ATMGES_ENSMEAN dimy | awk '{print $2}')} - LONB_ENKF=${LONB_ENKF:-$($NEMSIOGET $ATMGES_ENSMEAN dimx | awk '{print $2}')} - use_gfs_ncio=".false." - use_gfs_nemsio=".true." - paranc=${paranc:-".false."} -fi -LATA_ENKF=${LATA_ENKF:-$LATB_ENKF} -LONA_ENKF=${LONA_ENKF:-$LONB_ENKF} - -SATANGL=${SATANGL:-${FIXgsi}/global_satangbias.txt} -SATINFO=${SATINFO:-${FIXgsi}/global_satinfo.txt} -CONVINFO=${CONVINFO:-${FIXgsi}/global_convinfo.txt} -OZINFO=${OZINFO:-${FIXgsi}/global_ozinfo.txt} -SCANINFO=${SCANINFO:-${FIXgsi}/global_scaninfo.txt} -HYBENSINFO=${HYBENSINFO:-${FIXgsi}/global_hybens_info.l${LEVS_ENKF}.txt} -ANAVINFO=${ANAVINFO:-${FIXgsi}/global_anavinfo.l${LEVS_ENKF}.txt} -VLOCALEIG=${VLOCALEIG:-${FIXgsi}/vlocal_eig_l${LEVS_ENKF}.dat} - -ENKF_SUFFIX="s" -[[ $SMOOTH_ENKF = "NO" ]] && ENKF_SUFFIX="" - -################################################################################ -# Preprocessing -mkdata=NO -if [ ! -d $DATA ]; then - mkdata=YES - mkdir -p $DATA -fi -cd $DATA || exit 99 - -################################################################################ -# Fixed files -$NLN $SATANGL satbias_angle -$NLN $SATINFO satinfo -$NLN $SCANINFO scaninfo -$NLN $CONVINFO convinfo -$NLN $OZINFO ozinfo -$NLN $HYBENSINFO hybens_info -$NLN $ANAVINFO anavinfo -$NLN $VLOCALEIG vlocal_eig.dat - -# Bias correction coefficients based on the ensemble mean -$NLN $COMOUT_ANL_ENS/$GBIASe satbias_in - -################################################################################ - -if [ $USE_CFP = "YES" ]; then - [[ -f $DATA/untar.sh ]] && rm $DATA/untar.sh - [[ -f $DATA/mp_untar.sh ]] && rm $DATA/mp_untar.sh - set +x - cat > $DATA/untar.sh << EOFuntar -#!/bin/sh -memchar=\$1 -flist="$CNVSTAT $OZNSTAT $RADSTAT" -for ftype in \$flist; do - if [ \$memchar = "ensmean" ]; then - fname=$COMOUT_ANL_ENS/\${ftype}.ensmean - else - fname=$COMOUT_ANL_ENS/\$memchar/\$ftype - fi - tar -xvf \$fname -done -EOFuntar - set -x - chmod 755 $DATA/untar.sh -fi - -################################################################################ -# Ensemble guess, observational data and analyses/increments - -flist="$CNVSTAT $OZNSTAT $RADSTAT" -if [ $USE_CFP = "YES" ]; then - echo "$nm $DATA/untar.sh ensmean" | tee -a $DATA/mp_untar.sh - if [ ${CFP_MP:-"NO"} = "YES" ]; then - nm=$((nm+1)) - fi -else - for ftype in $flist; do - fname=$COMOUT_ANL_ENS/${ftype}.ensmean - tar -xvf $fname - done -fi -nfhrs=`echo $IAUFHRS_ENKF | sed 's/,/ /g'` -for imem in $(seq 1 $NMEM_ENKF); do - memchar="mem"$(printf %03i $imem) - if [ $lobsdiag_forenkf = ".false." ]; then - if [ $USE_CFP = "YES" ]; then - echo "$nm $DATA/untar.sh $memchar" | tee -a $DATA/mp_untar.sh - if [ ${CFP_MP:-"NO"} = "YES" ]; then - nm=$((nm+1)) - fi - else - for ftype in $flist; do - fname=$COMOUT_ANL_ENS/$memchar/$ftype - tar -xvf $fname - done - fi - fi - mkdir -p $COMOUT_ANL_ENSFSOI/$memchar - for FHR in $nfhrs; do - $NLN $COMIN_GES_ENS/$memchar/${GPREFIX}atmf00${FHR}${ENKF_SUFFIX}${GSUFFIX} sfg_${CDATE}_fhr0${FHR}_${memchar} - if [ $cnvw_option = ".true." ]; then - $NLN $COMIN_GES_ENS/$memchar/${GPREFIX}sfcf00${FHR}${GSUFFIX} sfgsfc_${CDATE}_fhr0${FHR}_${memchar} - fi - if [ $FHR -eq 6 ]; then - if [ $DO_CALC_INCREMENT = "YES" ]; then - #$NLN $COMOUT_ANL_ENS/$memchar/${APREFIX}atmanl${ASUFFIX} sanl_${CDATE}_fhr0${FHR}_${memchar} - $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atmanl${ASUFFIX} sanl_${CDATE}_fhr0${FHR}_ni${memchar} - else - #$NLN $COMOUT_ANL_ENS/$memchar/${APREFIX}atminc${ASUFFIX} incr_${CDATE}_fhr0${FHR}_${memchar} - $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atminc${ASUFFIX} incr_${CDATE}_fhr0${FHR}_ni${memchar} - fi - else - if [ $DO_CALC_INCREMENT = "YES" ]; then - #$NLN $COMOUT_ANL_ENS/$memchar/${APREFIX}atma00${FHR}${ASUFFIX} sanl_${CDATE}_fhr0${FHR}_${memchar} - $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atma00${FHR}${ASUFFIX} sanl_${CDATE}_fhr0${FHR}_ni${memchar} - else - #$NLN $COMOUT_ANL_ENS/$memchar/${APREFIX}atmi00${FHR}${ASUFFIX} incr_${CDATE}_fhr0${FHR}_${memchar} - $NLN $COMOUT_ANL_ENSFSOI/$memchar/${APREFIX}atmi00${FHR}${ASUFFIX} incr_${CDATE}_fhr0${FHR}_ni${memchar} - fi - fi - done -done - -# Ensemble mean guess -for FHR in $nfhrs; do - $NLN $COMIN_GES_ENS/${GPREFIX}atmf00${FHR}.ensmean${GSUFFIX} sfg_${CDATE}_fhr0${FHR}_ensmean - if [ $cnvw_option = ".true." ]; then - $NLN $COMIN_GES_ENS/${GPREFIX}sfcf00${FHR}.ensmean${GSUFFIX} sfgsfc_${CDATE}_fhr0${FHR}_ensmean - fi -done - -if [ $USE_CFP = "YES" ]; then - chmod 755 $DATA/mp_untar.sh - ncmd=$(cat $DATA/mp_untar.sh | wc -l) - if [ $ncmd -gt 0 ]; 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 - fi -fi - -# Keep osense data in comrot -$NLN $COMOUT_ANL_ENSFSOI/osense_${CDATE}_init.dat osense_${CDATE}.dat - -################################################################################ -# Create global_enkf namelist -cat > enkf.nml << EOFnml -&nam_enkf - datestring="$CDATE",datapath="$DATA/", - analpertwtnh=${analpertwt},analpertwtsh=${analpertwt},analpertwttr=${analpertwt}, - covinflatemax=1.e2,covinflatemin=1,pseudo_rh=.true.,iassim_order=0, - corrlengthnh=${corrlength},corrlengthsh=${corrlength},corrlengthtr=${corrlength}, - lnsigcutoffnh=${lnsigcutoff},lnsigcutoffsh=${lnsigcutoff},lnsigcutofftr=${lnsigcutoff}, - lnsigcutoffpsnh=${lnsigcutoff},lnsigcutoffpssh=${lnsigcutoff},lnsigcutoffpstr=${lnsigcutoff}, - lnsigcutoffsatnh=${lnsigcutoff},lnsigcutoffsatsh=${lnsigcutoff},lnsigcutoffsattr=${lnsigcutoff}, - obtimelnh=1.e30,obtimelsh=1.e30,obtimeltr=1.e30, - saterrfact=1.0,numiter=0, - sprd_tol=1.e30,paoverpb_thresh=0.98, - nlons=$LONA_ENKF,nlats=$LATA_ENKF,nlevs=$LEVS_ENKF,nanals=$NMEM_ENKF, - deterministic=.true.,sortinc=.true.,lupd_satbiasc=.false., - reducedgrid=${reducedgrid},readin_localization=${readin_localization_enkf}., - use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},imp_physics=$imp_physics,lupp=$lupp, - univaroz=.false.,adp_anglebc=.true.,angord=4,use_edges=.false.,emiss_bc=.true., - letkf_flag=${letkf_flag},nobsl_max=${nobsl_max},denkf=${denkf},getkf=${getkf}., - nhr_anal=${IAUFHRS_ENKF},nhr_state=${IAUFHRS_ENKF},use_qsatensmean=.true., - lobsdiag_forenkf=$lobsdiag_forenkf, - write_spread_diag=$write_spread_diag, - modelspace_vloc=$modelspace_vloc, - use_correlated_oberrs=${use_correlated_oberrs}, - netcdf_diag=$netcdf_diag,cnvw_option=$cnvw_option, - paranc=$paranc,write_fv3_incr=$write_fv3_incr, - efsoi_cycling=.true., - $WRITE_INCR_ZERO - $NAM_ENKF -/ -&satobs_enkf - sattypes_rad(1) = 'amsua_n15', dsis(1) = 'amsua_n15', - sattypes_rad(2) = 'amsua_n18', dsis(2) = 'amsua_n18', - sattypes_rad(3) = 'amsua_n19', dsis(3) = 'amsua_n19', - sattypes_rad(4) = 'amsub_n16', dsis(4) = 'amsub_n16', - sattypes_rad(5) = 'amsub_n17', dsis(5) = 'amsub_n17', - sattypes_rad(6) = 'amsua_aqua', dsis(6) = 'amsua_aqua', - sattypes_rad(7) = 'amsua_metop-a', dsis(7) = 'amsua_metop-a', - sattypes_rad(8) = 'airs_aqua', dsis(8) = 'airs_aqua', - sattypes_rad(9) = 'hirs3_n17', dsis(9) = 'hirs3_n17', - sattypes_rad(10)= 'hirs4_n19', dsis(10)= 'hirs4_n19', - sattypes_rad(11)= 'hirs4_metop-a', dsis(11)= 'hirs4_metop-a', - sattypes_rad(12)= 'mhs_n18', dsis(12)= 'mhs_n18', - sattypes_rad(13)= 'mhs_n19', dsis(13)= 'mhs_n19', - sattypes_rad(14)= 'mhs_metop-a', dsis(14)= 'mhs_metop-a', - sattypes_rad(15)= 'goes_img_g11', dsis(15)= 'imgr_g11', - sattypes_rad(16)= 'goes_img_g12', dsis(16)= 'imgr_g12', - sattypes_rad(17)= 'goes_img_g13', dsis(17)= 'imgr_g13', - sattypes_rad(18)= 'goes_img_g14', dsis(18)= 'imgr_g14', - sattypes_rad(19)= 'goes_img_g15', dsis(19)= 'imgr_g15', - sattypes_rad(20)= 'avhrr_n18', dsis(20)= 'avhrr3_n18', - sattypes_rad(21)= 'avhrr_metop-a', dsis(21)= 'avhrr3_metop-a', - sattypes_rad(22)= 'avhrr_n19', dsis(22)= 'avhrr3_n19', - sattypes_rad(23)= 'amsre_aqua', dsis(23)= 'amsre_aqua', - sattypes_rad(24)= 'ssmis_f16', dsis(24)= 'ssmis_f16', - sattypes_rad(25)= 'ssmis_f17', dsis(25)= 'ssmis_f17', - sattypes_rad(26)= 'ssmis_f18', dsis(26)= 'ssmis_f18', - sattypes_rad(27)= 'ssmis_f19', dsis(27)= 'ssmis_f19', - sattypes_rad(28)= 'ssmis_f20', dsis(28)= 'ssmis_f20', - sattypes_rad(29)= 'sndrd1_g11', dsis(29)= 'sndrD1_g11', - sattypes_rad(30)= 'sndrd2_g11', dsis(30)= 'sndrD2_g11', - sattypes_rad(31)= 'sndrd3_g11', dsis(31)= 'sndrD3_g11', - sattypes_rad(32)= 'sndrd4_g11', dsis(32)= 'sndrD4_g11', - sattypes_rad(33)= 'sndrd1_g12', dsis(33)= 'sndrD1_g12', - sattypes_rad(34)= 'sndrd2_g12', dsis(34)= 'sndrD2_g12', - sattypes_rad(35)= 'sndrd3_g12', dsis(35)= 'sndrD3_g12', - sattypes_rad(36)= 'sndrd4_g12', dsis(36)= 'sndrD4_g12', - sattypes_rad(37)= 'sndrd1_g13', dsis(37)= 'sndrD1_g13', - sattypes_rad(38)= 'sndrd2_g13', dsis(38)= 'sndrD2_g13', - sattypes_rad(39)= 'sndrd3_g13', dsis(39)= 'sndrD3_g13', - sattypes_rad(40)= 'sndrd4_g13', dsis(40)= 'sndrD4_g13', - sattypes_rad(41)= 'sndrd1_g14', dsis(41)= 'sndrD1_g14', - sattypes_rad(42)= 'sndrd2_g14', dsis(42)= 'sndrD2_g14', - sattypes_rad(43)= 'sndrd3_g14', dsis(43)= 'sndrD3_g14', - sattypes_rad(44)= 'sndrd4_g14', dsis(44)= 'sndrD4_g14', - sattypes_rad(45)= 'sndrd1_g15', dsis(45)= 'sndrD1_g15', - sattypes_rad(46)= 'sndrd2_g15', dsis(46)= 'sndrD2_g15', - sattypes_rad(47)= 'sndrd3_g15', dsis(47)= 'sndrD3_g15', - sattypes_rad(48)= 'sndrd4_g15', dsis(48)= 'sndrD4_g15', - sattypes_rad(49)= 'iasi_metop-a', dsis(49)= 'iasi_metop-a', - sattypes_rad(50)= 'seviri_m08', dsis(50)= 'seviri_m08', - sattypes_rad(51)= 'seviri_m09', dsis(51)= 'seviri_m09', - sattypes_rad(52)= 'seviri_m10', dsis(52)= 'seviri_m10', - sattypes_rad(53)= 'seviri_m11', dsis(53)= 'seviri_m11', - sattypes_rad(54)= 'amsua_metop-b', dsis(54)= 'amsua_metop-b', - sattypes_rad(55)= 'hirs4_metop-b', dsis(55)= 'hirs4_metop-b', - sattypes_rad(56)= 'mhs_metop-b', dsis(56)= 'mhs_metop-b', - sattypes_rad(57)= 'iasi_metop-b', dsis(57)= 'iasi_metop-b', - sattypes_rad(58)= 'avhrr_metop-b', dsis(58)= 'avhrr3_metop-b', - sattypes_rad(59)= 'atms_npp', dsis(59)= 'atms_npp', - sattypes_rad(60)= 'atms_n20', dsis(60)= 'atms_n20', - sattypes_rad(61)= 'cris_npp', dsis(61)= 'cris_npp', - sattypes_rad(62)= 'cris-fsr_npp', dsis(62)= 'cris-fsr_npp', - sattypes_rad(63)= 'cris-fsr_n20', dsis(63)= 'cris-fsr_n20', - sattypes_rad(64)= 'gmi_gpm', dsis(64)= 'gmi_gpm', - sattypes_rad(65)= 'saphir_meghat', dsis(65)= 'saphir_meghat', - sattypes_rad(66)= 'amsua_metop-c', dsis(66)= 'amsua_metop-c', - sattypes_rad(67)= 'mhs_metop-c', dsis(67)= 'mhs_metop-c', - sattypes_rad(68)= 'ahi_himawari8', dsis(68)= 'ahi_himawari8', - sattypes_rad(69)= 'abi_g16', dsis(69)= 'abi_g16', - sattypes_rad(70)= 'abi_g17', dsis(70)= 'abi_g17', - $SATOBS_ENKF -/ -&ozobs_enkf - sattypes_oz(1) = 'sbuv2_n16', - sattypes_oz(2) = 'sbuv2_n17', - sattypes_oz(3) = 'sbuv2_n18', - sattypes_oz(4) = 'sbuv2_n19', - sattypes_oz(5) = 'omi_aura', - sattypes_oz(6) = 'gome_metop-a', - sattypes_oz(7) = 'gome_metop-b', - sattypes_oz(8) = 'mls30_aura', - sattypes_oz(9) = 'ompsnp_npp', - sattypes_oz(10) = 'ompstc8_npp', - $OZOBS_ENKF -/ -EOFnml - -################################################################################ -# Run enkf update - -export OMP_NUM_THREADS=$NTHREADS_ENKF -export pgm=$ENKFEXEC -. prep_step - -$NCP $ENKFEXEC $DATA -$APRUN_ENKF ${DATA}/$(basename $ENKFEXEC) 1>stdout 2>stderr -rc=$? - -export ERR=$rc -export err=$ERR -$ERRSCRIPT || exit 2 - -# save for EFSOI task (still needed?) -$NCP $COMOUT_ANL_ENS/$GBIASe $COMOUT_ANL_ENSFSOI - -# Cat runtime output files. -cat stdout stderr > $COMOUT_ANL_ENSFSOI/$ENKFSTAT - -################################################################################ -# Postprocessing -######## AFE remove after testing -#cp -r $DATA /scratch1/NCEPDEV/stmp4/Andrew.Eichmann/efsoi -######## AFE -cd $pwd -[[ $mkdata = "YES" ]] && rm -rf $DATA -set +x -if [ $VERBOSE = "YES" ]; then - echo $(date) EXITING $0 with return code $err >&2 -fi -exit $err diff --git a/src/enkf/enkf.f90 b/src/enkf/enkf.f90 index 8ea45a5b67..248fe87717 100644 --- a/src/enkf/enkf.f90 +++ b/src/enkf/enkf.f90 @@ -129,9 +129,8 @@ module enkf iassim_order,sortinc,deterministic,numiter,nlevs,& zhuberleft,zhuberright,varqc,lupd_satbiasc,huber,univaroz,& covl_minfact,covl_efold,nbackgrounds,nhr_anal,fhr_assim,& - iseed_perturbed_obs,lupd_obspace_serial,efsoi_cycling,& + iseed_perturbed_obs,lupd_obspace_serial,fso_cycling,& neigv,vlocal_evecs,denkf - use radinfo, only: npred,nusis,nuchan,jpch_rad,predx use radbias, only: apply_biascorr, update_biascorr use gridinfo, only: nlevs_pres @@ -826,7 +825,7 @@ subroutine enkf_update() ! Gathering analysis perturbations ! in observation space for EFSO -if(efsoi_cycling) then +if(fso_cycling) then if(nproc /= 0) then call mpi_send(anal_obchunk,numobsperproc(nproc+1)*nanals,mpi_real,0, & 1,mpi_comm_world,ierr) diff --git a/src/enkf/enkf_main.f90 b/src/enkf/enkf_main.f90 index 1a38dd7525..be7df09608 100644 --- a/src/enkf/enkf_main.f90 +++ b/src/enkf/enkf_main.f90 @@ -76,7 +76,7 @@ program enkf_main ! reads namelist parameters. use params, only : read_namelist,cleanup_namelist,letkf_flag,readin_localization,lupd_satbiasc,& numiter, nanals, lupd_obspace_serial, write_spread_diag, & - lobsdiag_forenkf, netcdf_diag, efsoi_cycling, ntasks_io + lobsdiag_forenkf, netcdf_diag, fso_cycling, ntasks_io ! mpi functions and variables. use mpisetup, only: mpi_initialize, mpi_initialize_io, mpi_cleanup, nproc, & mpi_wtime @@ -183,7 +183,7 @@ program enkf_main ! Initialization for writing ! observation sensitivity files - if(efsoi_cycling) call init_ob_sens() + if(fso_cycling) call init_ob_sens() ! read in vertical profile of horizontal and vertical localization length ! scales, set values for each ob. @@ -216,7 +216,7 @@ program enkf_main ! Output non-inflated ! analyses for FSO - if(efsoi_cycling) then + if(fso_cycling) then no_inflate_flag=.true. t1 = mpi_wtime() call gather_chunks() @@ -240,7 +240,7 @@ program enkf_main endif ! print EFSO sensitivity i/o on root task. - if(efsoi_cycling) call print_ob_sens() + if(fso_cycling) call print_ob_sens() ! print innovation statistics for posterior on root task. if (nproc == 0 .and. numiter > 0) then @@ -268,7 +268,7 @@ program enkf_main call controlvec_cleanup() call loadbal_cleanup() - if(efsoi_cycling) call destroy_ob_sens() + if(fso_cycling) call destroy_ob_sens() call cleanup_namelist() ! write log file (which script can check to verify completion). diff --git a/src/enkf/enkf_obs_sensitivity.f90 b/src/enkf/enkf_obs_sensitivity.f90 index 66051c46b0..95aa71cf06 100644 --- a/src/enkf/enkf_obs_sensitivity.f90 +++ b/src/enkf/enkf_obs_sensitivity.f90 @@ -18,7 +18,6 @@ module enkf_obs_sensitivity ! destroy_ob_sens - Deallocate variables ! ! Variable Definitions: -! adloc_chunk - Coordinates of observation response ! obsense_kin - forecast sensitivity on each observations (kinetic energy) ! obsense_dry - forecast sensitivity on each observations (dry total energy) ! obsense_moist - forecast sensitivity on each observations (moist total energy) @@ -33,13 +32,13 @@ module enkf_obs_sensitivity use mpisetup, only: mpi_real4,mpi_sum,mpi_comm_io,mpi_in_place,numproc,nproc,& mpi_integer,mpi_wtime,mpi_status,mpi_real8,mpi_max,mpi_realkind use kinds, only: r_single,r_kind,r_double,i_kind -use params, only: efsoi_flag,latbound,nlevs,nanals,datestring, & +use params, only: fso_calculate,latbound,nlevs,nanals,datestring, & lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh, & lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh, & lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh, & corrlengthnh,corrlengthtr,corrlengthsh, & obtimelnh,obtimeltr,obtimelsh,letkf_flag, & - nbackgrounds,adrate,eft + nbackgrounds use constants, only: zero,one,half,rearth,pi,deg2rad,rad2deg use enkf_obsmod, only: nobstot,nobs_conv,nobs_oz,nobs_sat,obtype,obloclat, & obloclon,obpress,indxsat,oberrvar,stattype,obtime,ob, & @@ -49,8 +48,8 @@ module enkf_obs_sensitivity use convinfo, only: convinfo_read,init_convinfo use ozinfo, only: ozinfo_read,init_oz use radinfo, only: radinfo_read,jpch_rad,nusis,nuchan,npred -!use gridinfo_efsoi, only: latsgrd,lonsgrd,nlevs_pres,npts,id_u,id_v -use loadbal!, only: indxproc,grdloc_chunk,numptsperproc,npts_max,kdtree_grid +use gridinfo, only: latsgrd,lonsgrd,nlevs_pres,npts +use loadbal, only: indxproc,grdloc_chunk,numptsperproc,npts_max,kdtree_grid use covlocal, only: latval use kdtree2_module, only: kdtree2_create @@ -58,10 +57,9 @@ module enkf_obs_sensitivity private public init_ob_sens,destroy_ob_sens,print_ob_sens,read_ob_sens,& - obsense_kin,obsense_dry,obsense_moist,adloc_chunk + obsense_kin,obsense_dry,obsense_moist real(r_kind),allocatable,dimension(:) :: obsense_kin,obsense_dry,obsense_moist -real(r_single),allocatable,dimension(:,:) :: adloc_chunk ! Structure for observation sensitivity information output type obsense_header @@ -99,7 +97,6 @@ module enkf_obs_sensitivity contains - subroutine init_ob_sens !$$$ subprogram documentation block ! . . . . @@ -132,7 +129,6 @@ subroutine init_ob_sens end subroutine init_ob_sens - subroutine read_ob_sens !$$$ subprogram documentation block ! . . . . @@ -192,15 +188,8 @@ subroutine read_ob_sens write(6,*) 'READ_OBSENSE: number of members is not correct.',nanals,inhead%nanals call stop2(26) end if - + ! nobstot=nobsgood if(nproc == 0) write(6,*) 'total number of obs ',nobstot - if(nproc == 0) write(6,*) 'total number of conv obs ',nobs_conv - if(nproc == 0) write(6,*) 'total number of oz obs',nobs_oz - if(nproc == 0) write(6,*) 'total number of sat obs',nobs_sat - if(nproc == 0) write(6,*) 'npred=',inhead%npred - if(nproc == 0) write(6,*) 'idate=',inhead%idate - if(nproc == 0) write(6,*) 'nanals=',inhead%nanals - ! Allocate arrays allocate(obfit_prior(nobstot)) allocate(obsprd_prior(nobstot)) @@ -343,7 +332,7 @@ subroutine print_ob_sens type(obsense_info) :: outdata iunit = 10 ! Gather observation sensitivity informations to the root - if(efsoi_flag) then + if(fso_calculate) then allocate(recbuf(nobstot)) call mpi_reduce(obsense_kin,recbuf,nobstot,mpi_realkind,mpi_sum,0, & & mpi_comm_world,ierr) @@ -431,7 +420,7 @@ subroutine print_ob_sens outdata%stattype = stattype(nob) outdata%obtype = obtype(nob) outdata%indxsat = 0 - if(efsoi_flag) then + if(fso_calculate) then outdata%osense_kin = real(obsense_kin(nob),r_single) outdata%osense_dry = real(obsense_dry(nob),r_single) outdata%osense_moist = real(obsense_moist(nob),r_single) @@ -442,7 +431,7 @@ subroutine print_ob_sens end if tmpanal_ob(1:nanals) = real(anal_ob(1:nanals,nob),r_single) write(iunit) outdata,tmpanal_ob - if(.not. efsoi_flag) cycle + if(.not. fso_calculate) cycle ! Sum up nob_conv(iobtyp,ireg) = nob_conv(iobtyp,ireg) + 1 sumsense_conv(iobtyp,ireg,stkin) = sumsense_conv(iobtyp,ireg,stkin) & @@ -459,7 +448,7 @@ subroutine print_ob_sens rate_conv(iobtyp,ireg,stmoist) = rate_conv(iobtyp,ireg,stmoist) + one end do ! print out - if(efsoi_flag) then + if(fso_calculate) then print *,'observation impact for conventional obs' print *,'region, obtype, nobs, dJ, positive rate[%]:' do iobtyp=1,8 @@ -500,9 +489,9 @@ subroutine print_ob_sens outdata%oberrvar_orig = real(oberrvar_orig(nob),r_single) outdata%stattype = stattype(nob) outdata%obtype = obtype(nob) - outdata%indxsat = nchan + outdata%indxsat = nuchan(nchan) tmpbiaspreds(1:npred+1) = real(biaspreds(1:npred+1,nn),r_single) - if(efsoi_flag) then + if(fso_calculate) then outdata%osense_kin = real(obsense_kin(nob),r_single) outdata%osense_dry = real(obsense_dry(nob),r_single) outdata%osense_moist = real(obsense_moist(nob),r_single) @@ -513,7 +502,7 @@ subroutine print_ob_sens end if tmpanal_ob(1:nanals) = real(anal_ob(1:nanals,nob),r_single) write(iunit) outdata,tmpanal_ob,tmpbiaspreds - if(.not. efsoi_flag) cycle + if(.not. fso_calculate) cycle ! Sum up if (oberrvar(nob) < 1.e10_r_kind .and. nchan > 0) then nob_sat(nchan) = nob_sat(nchan) + 1 @@ -532,7 +521,7 @@ subroutine print_ob_sens end if end do ! print out - if(efsoi_flag) then + if(fso_calculate) then print *,'observation impact for satellite brightness temp' print *,'instrument, channel #, nobs, dJ, positive rate[%]:' do nchan=1,jpch_rad @@ -576,7 +565,6 @@ subroutine destroy_ob_sens if(allocated(obsense_kin)) deallocate(obsense_kin) if(allocated(obsense_dry)) deallocate(obsense_dry) if(allocated(obsense_moist)) deallocate(obsense_moist) - if(allocated(adloc_chunk)) deallocate(adloc_chunk) return end subroutine destroy_ob_sens end module enkf_obs_sensitivity diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index d5bd08a606..f245afdf3f 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -40,7 +40,7 @@ module gridio ! language: f95 ! !$$$ - use constants, only: zero,one,cp,fv,rd,tiny_r_kind,max_varname_length,t0c,r0_05,constants_initialized + use constants, only: zero,one,cp,fv,rd,tiny_r_kind,max_varname_length,t0c,r0_05 use params, only: nlons,nlats,nlevs,use_gfs_nemsio,pseudo_rh, & cliptracers,datapath,imp_physics,use_gfs_ncio,cnvw_option, & nanals diff --git a/src/enkf/loadbal.f90 b/src/enkf/loadbal.f90 index 595ef20277..60e7df8111 100644 --- a/src/enkf/loadbal.f90 +++ b/src/enkf/loadbal.f90 @@ -444,6 +444,8 @@ subroutine scatter_chunks end subroutine scatter_chunks + + subroutine gather_chunks ! gather chunks into grdin to write out the ensemble members use controlvec, only: ncdim, grdin diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 5bcfa062b9..89688eb4cc 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -30,8 +30,6 @@ module params ! modulated ensembles), nobsl_max (for ob selection ! in LETKF and dfs_sort ! (for using DFS in LETKF ob selection). -! 2018-11-15 groff - Added ancillary parameters -! for EFSOI calculations ! ! attributes: ! language: f95 @@ -57,13 +55,6 @@ module params integer(i_kind), public, parameter :: nsatmax_oz = 100 character(len=20), public, dimension(nsatmax_rad) ::sattypes_rad, dsis character(len=20), public, dimension(nsatmax_oz) ::sattypes_oz -! EFSOI file type identifiers -integer(i_kind), public, parameter :: read_ensmean_forecast = 0 -integer(i_kind), public, parameter :: read_analysis_mean = 1 -integer(i_kind), public, parameter :: read_member_forecasts = 2 -integer(i_kind), public, parameter :: read_verification = 3 -! Analysis impact specific file type identifier -integer(i_kind), public, parameter :: read_member_analyses = 2 ! forecast times for first-guess forecasts to be updated (in hours) integer,dimension(7),public :: nhr_anal = (/6,-1,-1,-1,-1,-1,-1/) integer,dimension(7),public :: nhr_state = (/6,-1,-1,-1,-1,-1,-1/) @@ -87,10 +78,6 @@ module params character(len=120),dimension(7),public :: incfileprefixes ! analysis date string (YYYYMMDDHH) character(len=10), public :: datestring -! Hour for datestring -integer(i_kind), public :: datehr, gdatehr -! analysis filename, needed for EFSOI calcs -character(len=100), public :: andataname ! filesystem path to input files (first-guess, GSI diagnostic files). character(len=500),public :: datapath ! if deterministic=.true., the deterministic square-root filter @@ -105,8 +92,6 @@ module params nanals_per_iotask, ntasks_io integer(i_kind),public, allocatable, dimension(:) :: nanal1,nanal2 integer(i_kind),public :: nsats_rad,nsats_oz,imp_physics -integer(i_kind),public :: eft -integer(i_kind),public :: tar_minlev,tar_maxlev ! random seed for perturbed obs (deterministic=.false.) ! if zero, system clock is used. Also used when ! iassim_order=1 (random shuffling of obs for serial assimilation). @@ -122,8 +107,6 @@ module params real(r_single),public :: analpertwtnh_rtpp,analpertwtsh_rtpp,analpertwttr_rtpp real(r_single),public :: paoverpb_thresh,latbound,delat,p5delat,delatinv real(r_single),public :: latboundpp,latboundpm,latboundmp,latboundmm -real(r_single),public :: wmoist,adrate -real(r_single),public :: tar_minlat,tar_maxlat,tar_minlon,tar_maxlon real(r_single),public :: covl_minfact, covl_efold real(r_single),public :: covinflatenh,covinflatesh,covinflatetr,lnsigcovinfcutoff @@ -188,11 +171,6 @@ module params logical,public :: nmm_restart = .true. logical,public :: nmmb = .false. logical,public :: letkf_flag = .false. - -! EFSOI ancillary flag to determine -! type of impact estimate/calculation -logical,public :: forecast_impact = .true. - ! use brute force search in LETKF instead of kdtree logical,public :: letkf_bruteforce_search=.false. @@ -204,12 +182,11 @@ module params logical,public :: dfs_sort = .false. ! if true generate additional input files -! required for EFSOI calculations -logical,public :: efsoi_cycling = .false. +! required for EFSO calculations +logical,public :: fso_cycling = .false. -! Ancillary flag, applied only for -! EFSOI calculation applications -logical,public :: efsoi_flag = .false. +! if true perform efso calculations +logical,public :: fso_calculate = .false. ! if true, use ensemble mean qsat in definition of ! normalized humidity analysis variable (instead of @@ -258,14 +235,11 @@ module params paoverpb_thresh,latbound,delat,pseudo_rh,numiter,biasvar,& lupd_satbiasc,cliptracers,simple_partition,adp_anglebc,angord,& newpc4pred,nmmb,nhr_anal,nhr_state, fhr_assim,nbackgrounds,nstatefields, & - save_inflation,nobsl_max,lobsdiag_forenkf,netcdf_diag,forecast_impact,& + save_inflation,nobsl_max,lobsdiag_forenkf,netcdf_diag,& letkf_flag,massbal_adjust,use_edges,emiss_bc,iseed_perturbed_obs,npefiles,& getkf,getkf_inflation,denkf,modelspace_vloc,dfs_sort,write_spread_diag,& covinflatenh,covinflatesh,covinflatetr,lnsigcovinfcutoff,letkf_bruteforce_search,& - efsoi_cycling,efsoi_flag,imp_physics,lupp,cnvw_option,use_correlated_oberrs,& - eft,wmoist,adrate,andataname,& - gdatehr,datehr,& - tar_minlat,tar_maxlat,tar_minlon,tar_maxlon,tar_minlev,tar_maxlev,& + fso_cycling,fso_calculate,imp_physics,lupp,cnvw_option,use_correlated_oberrs,& fv3_native, paranc, nccompress, write_fv3_incr,incvars_to_zero namelist /nam_wrf/arw,nmm,nmm_restart namelist /nam_fv3/fv3fixpath,nx_res,ny_res,ntiles,l_pres_add_saved @@ -283,10 +257,6 @@ subroutine read_namelist() ! defaults ! time (analysis time YYYYMMDDHH) datestring = "0000000000" ! if 0000000000 will not be used. -! default analysis hour -datehr = 00 -! Initial hour for background forecasts -gdatehr = 00 ! corrlength (length for horizontal localization in km) ! this corresponding GSI parameter is s_ens_h. ! corrlength is the distance at which the Gaspari-Cohn @@ -377,9 +347,6 @@ subroutine read_namelist() nlevs = 0 ! number of ensemble members nanals = 0 -! nvars is numer of 3d variables to update. -! for hydrostatic models, typically 5 (u,v,T,q,ozone). -nvars = 5 ! background error variance for rad bias coeffs (used in radbias.f90) ! default is (old) GSI value. ! if negative, bias coeff error variace is set to -biasvar/N, where @@ -388,22 +355,6 @@ subroutine read_namelist() ! analysis error variance from the previous cycle is used instead ! (same as in the GSI). biasvar = 0.1_r_single -! Evaluation FT for EFSOI -eft = 24 -! Weigt for moist total energy norm (0 when dry total energy) -! applied in EFSOI calculation -wmoist = 0.0_r_single -! Advection coefficient for localization function -adrate = 0.0_r_single -! Name of analysis file at EFSOI evaluation time -andataname='' -! Target area for observation impact computation -tar_minlat = -90.0_r_single -tar_maxlat = 90.0_r_single -tar_minlon = 0.0_r_single -tar_maxlon = 360.0_r_single -tar_minlev = 0 -tar_maxlev = 0 ! factor to multiply sat radiance errors. saterrfact = 1._r_single @@ -761,23 +712,6 @@ subroutine read_namelist() corrlengthtr = corrlengthtr * 1.e3_r_single/rearth corrlengthsh = corrlengthsh * 1.e3_r_single/rearth -if(efsoi_cycling) then - letkf_flag = .false. -end if - -! convert targe area boundary into radians -tar_minlat = tar_minlat * deg2rad -tar_maxlat = tar_maxlat * deg2rad -tar_minlon = tar_minlon * deg2rad -tar_maxlon = tar_maxlon * deg2rad - -! use default vertical levels -tar_maxlev = min(nlevs,tar_maxlev) -if(tar_minlev < 1 .or. tar_maxlev < 1 .or. tar_maxlev < tar_minlev) then - tar_minlev = 1 - tar_maxlev = nlevs -end if - ! this var is .false. until this routine is called. params_initialized = .true. diff --git a/src/gsi/constants.f90 b/src/gsi/constants.f90 index da2368d70f..caf4331532 100644 --- a/src/gsi/constants.f90 +++ b/src/gsi/constants.f90 @@ -84,10 +84,6 @@ module constants public :: soilmoistmin public :: stndrd_atmos_ps -! ------ EFSOI relevant parameters -------- ! - public :: tref, pref - public :: constants_initialized - ! Declare derived constants integer(i_kind):: huge_i_kind integer(i_kind), parameter :: max_varname_length=32 @@ -113,8 +109,6 @@ module constants real(r_kind),parameter:: ttp = 2.7316e+2_r_kind ! temperature at h2o triple point (K) real(r_kind),parameter:: jcal = 4.1855e+0_r_kind ! joules per calorie () real(r_kind),parameter:: stndrd_atmos_ps = 1013.25e2_r_kind ! 1976 US standard atmosphere ps (Pa) - real(r_kind),parameter:: tref = 2.8000e+2_r_kind ! reference T for total energy - real(r_kind),parameter:: pref = 1.0000e+5_r_kind ! reference P for total energy ! Numeric constants @@ -289,9 +283,6 @@ module constants integer(i_kind),parameter:: i_missing=-9999 integer(r_kind),parameter:: r_missing=-9999._r_kind -! Constants initialized - logical :: constants_initialized = .true. - contains subroutine init_constants_derived diff --git a/util/EFSOI_Utilities/fix/global_anavinfo_efsoi.l127.txt b/util/EFSOI_Utilities/fix/global_anavinfo_efsoi.l127.txt deleted file mode 100644 index 3ee946e74a..0000000000 --- a/util/EFSOI_Utilities/fix/global_anavinfo_efsoi.l127.txt +++ /dev/null @@ -1,120 +0,0 @@ -correlated_observations:: -! isis method kreq kmult type cov_file -iasi_metop-a 2 0.0 1.0 sea Rcov_iasiasea -iasi_metop-b 2 0.0 1.0 sea Rcov_iasibsea -iasi_metop-a 2 0.0 1.0 land Rcov_iasialand -iasi_metop-b 2 0.0 1.0 land Rcov_iasibland -cris-fsr_n20 2 0.0 1.0 sea Rcov_crisn20 -cris-fsr_npp 2 0.0 1.0 sea Rcov_crisnpp -:: - -met_guess:: -!var level crtm_use desc orig_name - ps 1 -1 surface_pressure ps - z 1 -1 geopotential_height phis - u 127 2 zonal_wind u - v 127 2 meridional_wind v - div 127 -1 zonal_wind div - vor 127 -1 meridional_wind vor - tv 127 2 virtual_temperature tv - q 127 2 specific_humidity sphu - oz 127 2 ozone ozone - cw 127 10 cloud_condensate cw - ql 127 12 cloud_liquid ql - qi 127 12 cloud_ice qi -:: - -state_derivatives:: -!var level src - ps 1 met_guess - u 127 met_guess - v 127 met_guess - tv 127 met_guess - q 127 met_guess - oz 127 met_guess - cw 127 met_guess - prse 128 met_guess -:: - -state_tendencies:: -!var levels source - u 127 met_guess - v 127 met_guess - tv 127 met_guess - q 127 met_guess - cw 127 met_guess - oz 127 met_guess - prse 128 met_guess -:: - -state_vector:: -!var level itracer source funcof - u 127 0 met_guess u - v 127 0 met_guess v - tv 127 0 met_guess tv - tsen 127 0 met_guess tv,q - q 127 1 met_guess q - oz 127 1 met_guess oz - ql 127 1 met_guess ql - qi 127 1 met_guess qi - prse 128 0 met_guess prse - ps 1 0 met_guess prse - sst 1 0 met_guess sst -:: - -state_vector_efsoi:: -!var level itracer source funcof - u 127 0 met_guess u - v 127 0 met_guess v - tv 127 0 met_guess tv -! tsen 127 0 met_guess tv,q - q 127 1 met_guess q -! oz 127 1 met_guess oz -! ql 127 1 met_guess ql -! qi 127 1 met_guess qi -! prse 128 0 met_guess prse - ps 1 0 met_guess prse -! sst 1 0 met_guess sst -:: - - -control_vector_enkf:: -!var level itracer as/tsfc_sdv an_amp0 source funcof - u 127 0 1.00 -1.0 state u,v - v 127 0 1.00 -1.0 state u,v - ps 1 0 1.20 -1.0 state prse -!pst 1 0 1.20 -1.0 state prse,u,v - tv 127 0 1.50 -1.0 state tv - q 127 1 1.50 -1.0 state q - oz 127 1 2.00 -1.0 state oz -!sst 1 0 1.00 -1.0 state sst -!cw 127 1 1.00 -1.0 state cw -!stl 1 0 3.00 -1.0 motley sst -!sti 1 0 3.00 -1.0 motley sst -:: - -control_vector:: -!var level itracer as/tsfc_sdv an_amp0 source funcof - sf 127 0 1.00 -1.0 state u,v - vp 127 0 1.00 -1.0 state u,v - ps 1 0 1.20 -1.0 state prse - t 127 0 1.50 -1.0 state tv - q 127 1 1.50 -1.0 state q - oz 127 1 2.00 -1.0 state oz - sst 1 0 1.00 -1.0 state sst - cw 127 1 1.00 -1.0 state cw - stl 1 0 3.00 -1.0 motley sst - sti 1 0 3.00 -1.0 motley sst -:: - -! Following table shows the use of all four prescribed trace gas data. -! To turn off any one of any combination of trace gas input, add "!" -! in the first column of that trace gas name. To use all default -! trace gas profiles, just delete the following seven lines. -chem_guess:: -!var level itracer crtm_use type orig_name -!ch4 64 1 2 n/a ch4 - co2 127 1 0 n/a co2 -!co 64 1 2 n/a co -!n2o 64 1 2 n/a n2o -:: diff --git a/util/EFSOI_Utilities/scripts/README_fv3gfs b/util/EFSOI_Utilities/scripts/README_fv3gfs deleted file mode 100644 index 0b4b870904..0000000000 --- a/util/EFSOI_Utilities/scripts/README_fv3gfs +++ /dev/null @@ -1,34 +0,0 @@ - -This file, README_fv3gfs, describes the files in GSI/util/EFSOI_Utilities/scripts -under the branch EXP-efso_fv3. These files are expected to be merged into master -shortly (as of 2020-07-30) - - -getefsiofiles.sh: obtains files necessary for gdaseupd task in global workflow, -the modified version of which is the first step in the EFSOI process - -copyefsoifiles.sh: copies files from experiment to EFSOI staging directory after -gdasepos has run - -display_osense.py: reads and plots sensitivity data for single cycle after reading -from osense file generated by EFSOI executable - -osense.py: contains routine to read osense file - -convdata_codes.csv: codes needed to organize convetional data sources, derived from -webpage (with the help of tools in pandas): -https://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_2.htm -May need to be periodically updated - -The following are left from Dave Groff's development work and are not longer viable, -but are retained for historical reference to aid current development. Once the -EFSOI utility is fully running they can be removed - -efso_advance_ensemble.sh -efso_cleanup.sh -efso_cleanup_wrapper.sh -efso_file_assembly.sh -efso_run.sh -efso_wrapper.sh -README_EFSOI - diff --git a/util/EFSOI_Utilities/scripts/convdata_codes.csv b/util/EFSOI_Utilities/scripts/convdata_codes.csv deleted file mode 100644 index 05c0c7be5c..0000000000 --- a/util/EFSOI_Utilities/scripts/convdata_codes.csv +++ /dev/null @@ -1,92 +0,0 @@ -,type,message,description,note,code -0,111,SYNDAT,"SYNTHETIC (BOGUS) TROPICAL CYCLONE STORM CENTER (generated in SYNDAT_SYNDATA) - q, Pstn",Not created (switch not set in SYNDAT_SYNDATA parm cards),111 -1,112,,"PSEUDO MEAN SEA-LEVEL PRESSURE AT TROPICAL CYCLONE STORM CENTER (generated in GSI, does not appear in pre-analysis PREPBUFR files) - Pstn",Pstn used by assimilation When implemented. may appear in post-analysis PREPBUFR file,112 -2,120,ADPUPA,"RAWINSONDE - Tv, q, Pstn, sst","Tv, q, Pstn used by assimilation sst monitored by assimilation by switch in convinfo text file read by GBL-GSI",120 -3,122,ADPUPA,"CLASS SOUNDING - Tv, q, Pstn",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Entire report neither monitored nor assimilated - not in convinfo text file read by GBL-GSI,122 -4,126,RASSDA,RASS [FROM NOAA PROFILER NETWORK (NPN) OR MULTI-AGENCY PROFILER (MAP) NETWORK] - Tv,Tv flagged for non-use by assimilation due to missing obs error Tv monitored by assimilation by switch in convinfo text file read by GBL-GSI NPN data are no longer produced after 13 September 2017. Multiu-Agency data are still available.,126 -5,130,AIRCFT,AIREP AND PIREP AIRCRAFT - Ts,Ts used by assimilation,130 -6,131 R,AIRCFT,"AMDAR AIRCRAFT - Ts, q (E-AMDAR only)",q (E-AMDAR only) not considered by assimilation Ts used by assimilation,131 -7,132,ADPUPA,"FLIGHT-LEVEL RECONNAISSANCE AND PROFILE DROPSONDE - Tv, q, Pstn","q (non-U.S. drops, NASA global Hawk drops, all levels), Pstn (all drops) [and surface level Tv (all drops), surface level q (all drops)] flagged for non-use by assimilation by switch in PREPOBS_PREPDATA parm cards Tv (all drops, above surface; reccos), q [U.S. (NOAA Gulf Stream and P-3, and USAF) drops, above surface; reccos] used by assimilation q (non-U.S. drops, NASA global Hawk drops, all levels) and Tv (all drops, at surface) monitored by assimilation due to their being flagged by switch in PREPOBS_PREPDATA parm cards",132 -8,133 R,AIRCAR,"MDCRS ACARS AIRCRAFT - Ts, q","Ts, q used by assimilation",133 -9,134 R,AIRCFT,"TAMDAR AIRCRAFT - Ts, q","Ts, q monitored by assimilation by switch in convinfo text file read by GBL-GSI",134 -10,135 R,AIRCFT,CANADIAN AMDAR AIRCRAFT - Ts,Ts monitored by assimilation by switch in convinfo text file read by GBL-GSI,135 -11,150,SPSSMI,SSM/I SUPEROBED (1 DEGREE LAT/LON) FNMOC (OPERATIONAL) RAIN RATE (DMSP) - rr,rr monitored by assimilation by switch in pcpinfo text file read by GBL-GSI Currently the assimilation obtains this directly from the spssmi dump file rather than from PREPBUFR file. The SSM/I F-13 satellite went bad in November 2009 resulting in no data being processed upstream.  The empty dumps were turned off and the PREPOBS_PREPDATA parm cards were set to no longer process these data in October 2010.,150 -12,151,GOESND,"NESDIS 1x1 F-O-V CLOUD TOP PRESSURE, TEMPERATURE; CLOUD AMOUNT (GOES)",Not dumped,151 -13,152,SPSSMI,SSM/I SUPEROBED (1 DEGREE LAT/LON) NEURAL NET-3 PRECIPITABLE WATER OVER OCEAN (DMSP) - PWt,PWt flagged for non-use by assimilation due to missing obs error PWt monitored by assimilation by switch in convinfo text file read by GBL-GSI The SSM/I F-13 satellite went bad in November 2009 resulting in no data being processed upstream.  The empty dumps were turned off and the PREPOBS_PREPDATA parm cards were set to no longer process these data in October 2010.,152 -14,153,GPSIPW,GPS-INTEGRATED PRECIPITABLE WATER (GPS-IPW) - PWt,Currently only reports from ENI (mainly over U.S.) are encoded into PREPBUFR file. European GNSS reports are skipped. PWt flagged for non-use by assimilation due to missing obs error PWt monitored by assimilation by switch in convinfo text file read by GBL-GSI,153 -15,154,,reserved for GOES IMAGER SKY-COVER DATA used only in RTMA/URMA,,154 -16,156,GOESND,NESDIS 1x1 F-O-V 4-LAYER PRECIPITABLE WATER OVER LAND - CLEAR (GOES) - PWl,Not dumped,156 -17,157,GOESND,NESDIS 1x1 F-O-V 4-LAYER PRECIPITABLE WATER OVER LAND - CLOUDY (GOES) - PWl,Not dumped,157 -18,158,GOESND,NESDIS 1x1 F-O-V 4-LAYER PRECIPITABLE WATER OVER OCEAN - CLEAR (GOES) - PWl,Not dumped,158 -19,159,GOESND,NESDIS 1x1 F-O-V 4-LAYER PRECIPITABLE WATER OVER OCEAN - CLOUDY (GOES) - PWl,Not dumped,159 -20,164,GOESND,NESDIS 1x1 F-O-V RADIANCES OVER LAND - CLEAR (GOES) - Tb,Not dumped,164 -21,165,GOESND,NESDIS 1x1 F-O-V RADIANCES OVER LAND - CLOUDY (GOES) - Tb,Not dumped,165 -22,170,,"NACELLE - Tb, q",Not dumped,170 -23,171,,"TALL TOWER - Tb, q",Not dumped,171 -24,174,GOESND,NESDIS 1x1 F-O-V RADIANCES OVER OCEAN - CLEAR (GOES) - Tb,Not dumped,174 -25,175,GOESND,NESDIS 1x1 F-O-V RADIANCES OVER OCEAN - CLOUDY (GOES) - Tb,Not dumped,175 -26,180 (R - U.S. & JMA SHIPS),SFCSHP,"SURFACE MARINE WITH REPORTED STATION PRESSURE (SHIP, BUOY, C-MAN, TIDE GAUGE) - Tv, q, Pstn, sst","Tv, q, Pstn used by assimilation sst monitored by assimilation by switch in convinfo text file read by GBL-GSI",180 -27,181 (R - WMO Res 40 SYNOPS),ADPSFC,"SURFACE LAND [SYNOPTIC (fixed and mobile), METAR] WITH REPORTED STATION PRESSURE - Tv, q, Pstn, sst","Pstn used by assimilation Tv, q flagged for non-use by assimilation due to missing obs error Tv, q, sst monitored by assimilation by switch in convinfo text file read by GBL-GSI",181 -28,182,SFCSHP,"SPLASH-LEVEL DROPSONDE OVER OCEAN - Tv, q, Pstn","Tv, q, Pstn used by assimilation",182 -29,"183 (R - WMO Res 40 SYNOPS, U.S. & JMA SHIPS)","ADPSFC, SFCSHP","SURFACE MARINE (SHIP, BUOY, C-MAN, TIDE GAUGE) OR LAND [SYNOPTIC (fixed and mobile), METAR] WITH MISSING STATION PRESSURE - Tv, q, Pstn, sst","Tv, q, Pstn (entire report) flagged for non-use by assimilation due to missing obs error Tv, q, Pstn, sst monitored by assimilation by switch in convinfo text file read by GBL-GSI Altimeter setting is also missing (always the case for synoptic). Station pressure calculated from reported mean sea-level pressure and elevation via U.S. Standard Atmosphere approximation. Elevation is greater than 7.5 meters (if less than 7.5 meters, station pressure set equal to sea-level pressure and report type set to 181).",183 -30,187,ADPSFC,"SURFACE LAND (METAR) WITH MISSING STATION PRESSURE - Tv, q, Pstn, sst","Pstn used by assimilation Tv, q flagged for non-use by assimilation due to missing obs error Tv, q, sst monitored by assimilation by switch in convinfo text file read by GBL-GSI Altimeter setting is reported (never the case for synoptic). Station pressure calculated from reported altimeter setting and elevation.",187 -31,188 R,MSONET,"SURFACE MESONET - Tv, q, Pstn",Not dumped,188 -32,191,SFCBOG,AUSTRALIAN PAOB MEAN SEA-LEVEL PRESSURE BOGUS OVER OCEAN - Pstn,sfcbog dump file not read by PREPOBS_PREPDATA This data no longer produced after 17 August 2010,191 -33,192,ADPSFC,"SURFACE LAND SYNOPIC (fixed and mobile) WITH MISSING STATION PRESSURE AND MISSING SEA-LEVEL PRESSURE - Tv, q, Pstn",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Station pressure estimated from U.S. Standard Atmosphere approximation Pmsl and reported temperature and elevation.,192 -34,193,ADPSFC,"SURFACE LAND METAR WITH MISSING STATION PRESSURE, MISSING SEA-LEVEL PRESSURE AND MISSING ALTIMETER SETTING - Tv, q, Pstn",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Station pressure estimated from U.S. Standard Atmosphere approximation Pmsl and reported temperature and elevation.,193 -35,194,SFCSHP,"SURFACE MARINE (SHIP, BUOY, C-MAN, TIDE GAUGE) OR LAND (SYNOPTIC, METAR) WITH MISSING STATION PRESSURE AND MISSING SEA-LEVEL PRESSURE - Tv, q, Pstn",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Station pressure estimated from U.S. Standard Atmosphere approximation Pmsl and reported temperature and elevation.,194 -36,195,MSONET,"SURFACE MESONET WITH MISSING STATION PRESSURE AND MISSING ALTIMETER SETTING (SEA-LEVEL PRESSURE IS ALWAYS MISSING) - Tv, q, Pstn",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Station pressure estimated from U.S. Standard Atmosphere approximation Pmsl and reported temperature and elevation.,195 -37,210,SYNDAT,"SYNTHETIC (BOGUS) TROPICAL CYCLONE - u, v","u, v monitored by assimilation by switch in convinfo text file read by GBL-GSI Added to PREPBUFR file by later program SYNDAT_SYNDATA.",210 -38,220,ADPUPA,"RAWINSONDE - u, v (all levels), z (winds-by-height levels)","u, v used by assimilation p (vertical coordinate) calculated from z on winds-by-height levels.",220 -39,221,ADPUPA,"PIBAL - u,v,z","u, v used by assimilation p (vertical coordinate) calculated from z.",221 -40,222,ADPUPA,"CLASS SOUNDING - u, v",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Entire report neither monitored nor assimilated - not in convinfo text file read by GBL-GSI,222 -41,223,PROFLR,"NOAA PROFILER NETWORK (NPN) WIND PROFILER - u, v, z","u, v used by assimilation p (vertical coordinate) calculated from z. These data are no longer produced after 13 September 2017.",223 -42,224,VADWND,"NEXRAD VERTICAL AZIMUTH DISPLAY (VAD) from Radar Coded Message (subtype 1) - u, v, z","u, v used by assimilation p (vertical coordinate) calculated from z.",224 -43,227,PROFLR,"MULTI-AGENCY PROFILER (MAP) AND ACOUSTIC SOUNDER (SODAR) - u, v, z",Entire report tossed by switch in PREPOBS_PREPDATA parm cards,227 -44,228,PROFLR,"JAPANESE METEOROLOGICAL AGENCY (JMA) WIND PROFILER - u, v, z","u, v (entire report) flagged for non-use by assimilation due to missing obs error u, v monitored by assimilation by switch in convinfo text file read by GBL-GSI",228 -45,229,PROFLR,"WIND PROFILER DECODED FROM PILOT (PIBAL) BULLETINS - u, v, z","u, v used by assimilation p (vertical coordinate) calculated from z.",229 -46,230,AIRCFT,"AIREP AND PIREP AIRCRAFT - u, v","u, v used by assimilation",230 -47,231 R,AIRCFT,"AMDAR AIRCRAFT - u, v","u, v used by assimilation",231 -48,232,ADPUPA,"FLIGHT-LEVEL RECONNAISSANCE AND PROFILE DROPSONDE - u, v","u, v used by assimilation",232 -49,233 R,AIRCAR,"MDCRS ACARS AIRCRAFT - u, v","u, v used by assimilation",233 -50,234 R,AIRCFT,"TAMDAR AIRCRAFT - u, v","u, v monitored by assimilation by switch in convinfo text file read by GBL-GSI",234 -51,235 R,AIRCFT,"CANADIAN AMDAR AIRCRAFT - u, v","u, v monitored by assimilation by switch in convinfo text file read by GBL-GSI",235 -52,240,SATWND,"NESDIS IR (SHORT-WAVE) CLOUD DRIFT (ALL LEVELS) (GOES) - u, v",The assimilation will obtain this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This will not be written into the PREPBUFR file.,240 -53,241,SATWND,"INDIA IR (LONG-WAVE) AND VISIBLE CLOUD DRIFT (ALL LEVELS) (INSAT, KALPANA) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",241 -54,242,SATWND,"JMA IR (LONG-WAVE) AND VISIBLE CLOUD DRIFT AT LEVELS BELOW 850 MB (GMS, MTSAT, HIMAWARI) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI. The GSI redefines report type 242 as JMA visible cloud drift at all levels.",242 -55,243,SATWND,"EUMETSAT IR (LONG-WAVE) AND VISIBLE CLOUD DRIFT AT LEVELS BELOW 850 MB (METEOSAT) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI. The GSI redefines report type 243 as EUMETSAT visible cloud drift at all levels.",243 -56,244,SATWND,"AVHRR/POES IR (LONG-WAVE) CLOUD DRIFT (ALL LEVELS) (NOAA, METOP) - u,v",The assimilation will obtain this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This will not be written into the PREPBUFR file.,244 -57,245,SATWND,"NESDIS IR (LONG-WAVE) CLOUD DRIFT (ALL LEVELS) (GOES) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",245 -58,246,SATWND,"NESDIS IMAGER WATER VAPOR (ALL LEVELS) - CLOUD TOP (GOES) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",246 -59,247,SATWND,"NESDIS IMAGER WATER VAPOR (ALL LEVELS) - DEEP LAYER (GOES) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file (it was never in the PREPBUFR file) (see Table 18).  Since it has both a missing obs error and is set to monitor in the convinfo file, it is now monitored by the GBL-GSI.",247 -60,248,SATWND,"NESDIS SOUNDER WATER VAPOR (ALL LEVELS) - CLOUD TOP (GOES) - u, v","If ever processed, the assimilation will obtain this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This will never be written into the PREPBUFR file.",248 -61,249,SATWND,"NESDIS SOUNDER WATER VAPOR (ALL LEVELS) - DEEP LAYER (GOES) - u, v","If ever processed, the assimilation will obtain this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This will never be written into the PREPBUFR file.",249 -62,250,SATWND,"JMA IMAGER WATER VAPOR (ALL LEVELS) - CLOUD TOP & DEEP LAYER (GMS, MTSAT, HIMAWARI) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",250 -63,251,SATWND,"NESDIS VISIBLE CLOUD DRIFT (ALL LEVELS) (GOES) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",251 -64,252,SATWND,"JMA IR (LONG-WAVE) AND VISIBLE CLOUD DRIFT AT LEVELS ABOVE 850 MB (GMS, MTSAT, HIMAWARI) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI. The GSI redefines report type 252 as JMA IR cloud drift at all levels.",252 -65,253,SATWND,"EUMETSAT IR (LONG-WAVE) AND VISIBLE CLOUD DRIFT AT LEVELS ABOVE 850 MB (METEOSAT) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI. The GSI redefines report type 253 as EUMETSAT IR cloud drift at all levels.",253 -66,254,SATWND,"EUMETSAT IMAGER WATER VAPOR (ALL LEVELS) - CLOUD TOP & DEEP LAYER (METEOSAT) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",254 -67,255,SATWND,NESDIS PICTURE TRIPLET CLOUD DRIFT (LOW LEVELS) (GOES),No longer produced by NESDIS,255 -68,256,SATWND,"INDIA IMAGER WATER VAPOR (ALL LEVELS) (INSAT, KALPANA) - u, v","If ever processed, the assimilation will obtain this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This will never be written into the PREPBUFR file.",256 -69,257,SATWND,"MODIS/POES IR (LONG-WAVE) CLOUD DRIFT (ALL LEVELS) (AQUA, TERRA) - u,v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",257 -70,258,SATWND,"MODIS/POES IMAGER WATER VAPOR (ALL LEVELS) - CLOUD TOP (AQUA, TERRA) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",258 -71,259,SATWND,"MODIS/POES IMAGER WATER VAPOR (ALL LEVELS) - DEEP LAYER (AQUA, TERRA) - u, v","Effective 5/22/2012, the assimilation obtains this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This is still written into the PREPBUFR file but is ignored by the GBL-GSI.",259 -72,260,SATWND,"VIIRS/POES IR (LONG-WAVE) CLOUD DRIFT (ALL LEVELS) (NPP) - u,v",The assimilation will obtain this directly from the satwnd dump file rather than from the PREPBUFR file (see Table 18).  This will not be written into the PREPBUFR file.,260 -73,270,,"NACELLE - u,v",Not dumped,270 -74,271,,"TALL TOWER -u,v",Not dumped,271 -75,280 (R - U.S. & JMA SHIPS),SFCSHP,"SURFACE MARINE WITH REPORTED STATION PRESSURE (SHIP, BUOY, C-MAN, TIDE GAUGE) - u, v","u, v used by assimilation",280 -76,281 (R - WMO Res 40 SYNOPS),ADPSFC,"SURFACE LAND [SYNOPTIC (fixed and mobile), METAR] WITH REPORTED STATION PRESSURE - u, v","u, v (entire report) flagged for non-use by assimilation due to missing obs error u, v monitored by assimilation by switch in convinfo text file read by GBL-GSI",281 -77,282,SFCSHP,"ATLAS BUOY - u, v (see % below)","u, v used by assimilation Reported station pressure and mean sea-level pressure BOTH missing. Station pressure is set to 1013 mb. Elevation is less than or equal to 7.5 meters.",282 -78,283,SPSSMI,"SSM/I SUPEROBED (1 DEGREE LAT/LON) NEURAL NET 3 WIND SPEED OVER OCEAN - u, v","u, v monitored by assimilation by switch in convinfo text file read by GBL-GSI Only wspd available so direction initially set to ZERO in PREPBUFR file, direction calculated from analysis and u, v recomputed in program PREPOBS_OIQCBUFR. Reported station pressure and mean sea-level pressure BOTH missing. Station pressure is set to 1013 mb. Elevation is less than or equal to 7.5 meters. The SSM/I F-13 satellite went bad in November 2009 resulting in no data being processed upstream.  The empty dumps were turned off and the PREPOBS_PREPDATA parm cards were set to no longer process these data in October 2010.",283 -79,"284 (R - WMO Res 40 SYNOPS, U.S. & JMA SHIPS)","ADPSFC, SFCSHP","SURFACE MARINE (SHIP, BUOY, C-MAN, TIDE GAUGE) OR LAND [SYNOPTIC (fixed and mobile), METAR] WITH MISSING STATION PRESSURE - u, v","u, v (entire report) flagged for non-use by assimilation due to missing obs error u, v monitored by assimilation by switch in convinfo text file read by GBL-GSI Altimeter setting is also missing (always the case for synoptic). Station pressure calculated from reported mean sea-level pressure and elevation via U.S. Standard Atmosphere approximation. Elevation is greater than 7.5 meters (if less than 7.5 meters, station pressure set equal to sea-level pressure and report type set to 281).",284 -80,285,QKSWND,"SUPEROBED (0.5 DEGREE LAT/LON) SCATTEROMETER WINDS OVER OCEAN (QUIKSCAT) - u, v",No longer produced,285 -81,286,ERS1DA,"SCATTEROMETER WINDS OVER OCEAN (ERS) - u, v",No longer produced,286 -82,287,ADPSFC,"SURFACE LAND (METAR) WITH MISSING STATION PRESSURE - u, v","u, v (entire report) flagged for non-use by assimilation due to missing obs error u, v monitored by assimilation by switch in convinfo text file read by GBL-GSI Altimeter setting is reported (never the case for synoptic). Station pressure calculated from reported altimeter setting and elevation.",287 -83,288 R,MSONET,"SURFACE MESONET - u, v",Not dumped,288 -84,289,WDSATR,"SUPEROBED (1.0 DEGREE LAT/LON) SCATTEROMETER WINDS OVER OCEAN (WINDSAT) - u,v","u, v used by assimilation (currently not available due to format change in raw files)",289 -85,290,ASCATW,"NON-SUPEROBED SCATTEROMETER WINDS OVER OCEAN (ASCAT) - u,v (50 km resolution)","METOP-2(A) u, v used by assimilation METOP-1(B) u,v monitored by assimilation by switch in convinfo text file read by GBL-GSI",290 -86,291,,"NON-SUPEROBED SCATTEROMETER WINDS OVER OCEAN (OSCAT) - u,v","When available, the assimilation will obtain this directly from the oscatw dump file rather than from the PREPBUFR file (see Table 18).  This will not be written into the PREPBUFR file. This type will never be available. Instrument failed on 20 February 2014.",291 -87,292,ADPSFC,"SURFACE LAND SYNOPIC (fixed and mobile) WITH MISSING STATION PRESSURE AND MISSING SEA-LEVEL PRESSURE - u,v",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Station pressure estimated from U.S. Standard Atmosphere approximation Pmsl and reported temperature and elevation.,292 -88,293,ADPSFC,"SURFACE LAND METAR WITH MISSING STATION PRESSURE, MISSING SEA-LEVEL PRESSURE AND MISSING ALTIMETER SETTING - u,v",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Station pressure estimated from U.S. Standard Atmosphere approximation Pmsl and reported temperature and elevation.,293 -89,294,SFCSHP,"SURFACE MARINE (SHIP, BUOY, C-MAN, TIDE GAUGE) OR LAND (SYNOPTIC, METAR) WITH MISSING STATION PRESSURE AND MISSING SEA-LEVEL PRESSURE - u,v",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Station pressure estimated from U.S. Standard Atmosphere approximation Pmsl and reported temperature and elevation.,294 -90,295,MSONET,"SURFACE MESONET WITH MISSING STATION PRESSURE AND MISSING ALTIMETER SETTING (SEA-LEVEL PRESSURE IS ALWAYS MISSING) - u,v",Entire report tossed by switch in PREPOBS_PREPDATA parm cards Station pressure estimated from U.S. Standard Atmosphere approximation Pmsl and reported temperature and elevation.,295 diff --git a/util/EFSOI_Utilities/scripts/copyefsoifiles.sh b/util/EFSOI_Utilities/scripts/copyefsoifiles.sh deleted file mode 100755 index 7dc5bb6da6..0000000000 --- a/util/EFSOI_Utilities/scripts/copyefsoifiles.sh +++ /dev/null @@ -1,107 +0,0 @@ -#!/bin/sh -x - - -EXPDIR=/scratch1/NCEPDEV/stmp4/Andrew.Eichmann/nov2019b -STGDIR=/scratch1/NCEPDEV/stmp4/Andrew.Eichmann/efsoitest2/staging -#STGDIR=/scratch1/NCEPDEV/stmp4/Andrew.Eichmann/tmp/ -#EDATE=2019111800 -EDATE=2019111818 -NENS=20 -NDATE=/scratch1/NCEPDEV/global/Fanglin.Yang/save/VRFY/vsdb/nwprod/util/exec/ndate - -PDATE=`${NDATE} -6 ${EDATE}` -VDATE=`${NDATE} +24 ${EDATE}` - -EDATECYC="${EDATE:8:2}" -PDATECYC="${PDATE:8:2}" -VDATECYC="${VDATE:8:2}" - -EDATEDIR=enkfgdas."${EDATE:0:8}" -PDATEDIR=enkfgdas."${PDATE:0:8}" -VDATEDIR=enkfgdas."${VDATE:0:8}" - - -echo $PDATE -echo $VDATE -echo $EDATECYC -echo $EDATEDIR - -# set up 30 hour forecast - -cd $STGDIR - -if [ -d "$PDATEDIR" ]; then - echo "$PDATEDIR exists." -else - echo "$PDATEDIR does not exist, creating it" - mkdir $PDATEDIR -fi - -cd $PDATEDIR - -if [ -d "$PDATECYC" ]; then - echo "$PDATECYC exists." -else - echo "$PDATECYC does not exist, creating it" - mkdir $PDATECYC -fi - -cd $PDATECYC - -cp -vi $EXPDIR/$PDATEDIR/$PDATECYC/gdas.t${PDATECYC}z.atmf030.ensmean.* . -cp -vi $EXPDIR/$PDATEDIR/$PDATECYC/gdas.t${PDATECYC}z.atmf006.ensmean.* . - - -# set up 24 hour forecast - -cd $STGDIR - -if [ -d "$EDATEDIR" ]; then - echo "$EDATEDIR exists." -else - echo "$EDATEDIR does not exist, creating it" - mkdir $EDATEDIR -fi - -cd $EDATEDIR - -if [ -d "$EDATECYC" ]; then - echo "$EDATECYC exists." -else - echo "$EDATECYC does not exist, creating it" - mkdir $EDATECYC -fi - -cd $EDATECYC - -cp -vi $EXPDIR/$EDATEDIR/$EDATECYC/gdas.t${EDATECYC}z.atmf024.ensmean.* . -cp -vi $EXPDIR/$EDATEDIR/$EDATECYC/gdas.t${EDATECYC}z.atmf006.ensmean.* . -#cp -vi $EXPDIR/$EDATEDIR/$EDATECYC/gdas.t${EDATECYC}z.atmanl.ensmean.* . -cp -vi $EXPDIR/$EDATEDIR/$EDATECYC/gdas.t${EDATECYC}z.abias_int.ensmean . - - - -imem=1 -while [[ $imem -le $NENS ]]; do - member="mem"`printf %03i $imem` - - if [ -d "$member" ]; then - echo "$member exists." - else - echo "$member does not exist, creating it" - mkdir $member - fi - - cd $member - - cp -vi $EXPDIR/$EDATEDIR/$EDATECYC/$member/gdas.t${EDATECYC}z.atmf024.nemsio . - cp -vi $EXPDIR/$EDATEDIR/$EDATECYC/$member/gdas.t${EDATECYC}z.atmf024s.nemsio . - - cd .. - - (( imem = $imem + 1 )) -done - - - - diff --git a/util/EFSOI_Utilities/scripts/display_osense.py b/util/EFSOI_Utilities/scripts/display_osense.py deleted file mode 100644 index ae6f6b012c..0000000000 --- a/util/EFSOI_Utilities/scripts/display_osense.py +++ /dev/null @@ -1,49 +0,0 @@ -import osense -import pandas as pd - - -( convdata, satdata )= osense.read_osense('osense_2017091000.dat.out') - -# creates a DataFrame with the mean of each satellite instrument -#meanbyobtype = satdata.groupby('obtype').mean() - -meanbyobtype=satdata[['obtype','osense_kin','osense_dry','osense_moist']].groupby('obtype').mean() -#cmeanbyobtype=convdata[['stattype','osense_kin','osense_dry','osense_moist']].groupby('stattype').mean() - -convcodes = pd.read_csv('convdata_codes.csv') - -# associate each data point with its source, by code -# it would be more efficient to take the mean of the observation sensitivities by -# code/stattype, but this way the mean is by the message column in the codes -# (ADPUPA, AIRCRAFT, etc) to consolidate for simpler graphing. Taking the mean -# by code/stattype would break it down by data source more -convbycodes=pd.merge(convdata,convcodes,how='left',left_on='stattype', right_on='code') - -convmean=convbycodes[['message','osense_kin','osense_dry','osense_moist']].groupby('message').mean() -alltheobtypes=pd.concat([meanbyobtype,convmean]) - -figuresize = (10,6) - -o_dry=alltheobtypes['osense_dry'].sort_values(ascending=False) -dry_plot=o_dry.plot.barh(title = '2017091000 mean obs sensitivity, dry', figsize=figuresize); -fig=dry_plot.get_figure() -fig.savefig('obsense_dry.png',bbox_inches='tight') - -fig.clear() - -o_moist=alltheobtypes['osense_moist'].sort_values(ascending=False) -moist_plot=o_moist.plot.barh(title = '2017091000 mean obs sensitivity, moist', figsize=figuresize); -fig=moist_plot.get_figure() -fig.savefig('obsense_moist.png',bbox_inches='tight') - -fig.clear() - -o_kin=alltheobtypes['osense_kin'].sort_values(ascending=False) -kin_plot=o_kin.plot.barh(title = '2017091000 mean obs sensitivity, kinetic', figsize=figuresize); -fig=kin_plot.get_figure() -fig.savefig('obsense_kin.png',bbox_inches='tight') - -fig.clear() - - - diff --git a/util/EFSOI_Utilities/scripts/getefsiofiles.sh b/util/EFSOI_Utilities/scripts/getefsiofiles.sh deleted file mode 100755 index 3c41b7722a..0000000000 --- a/util/EFSOI_Utilities/scripts/getefsiofiles.sh +++ /dev/null @@ -1,77 +0,0 @@ -#!/bin/sh -x - - -# Author: Andrew Eichmann -# Purpose: to obtain from the HPSS archive the minimal set of files needed to -# run the gdaseupd task of GFS workflow, intended for EFSOI processing -# Script assumes that these files have been archived during a run; this will -# probably require modifications to the archiving script in workflow. - -# This will have to be changed to handle nc4 files - -EXP=nov2019c -#EDATE=2019111800 -#EDATE=2019111818 -EDATE=2019111712 -NMEM_ENKF=20 # total number of ensemble members -#NMEM_ENKF=0 -NMEM_EARCGRP=10 # number of ensemble members in archiving group -NDATE=/scratch1/NCEPDEV/global/Fanglin.Yang/save/VRFY/vsdb/nwprod/util/exec/ndate - -PDATE=`${NDATE} -6 ${EDATE}` -VDATE=`${NDATE} +24 ${EDATE}` - -EDATECYC="${EDATE:8:2}" -PDATECYC="${PDATE:8:2}" -VDATECYC="${VDATE:8:2}" - -EDATEDIR=enkfgdas."${EDATE:0:8}" -PDATEDIR=enkfgdas."${PDATE:0:8}" -VDATEDIR=enkfgdas."${VDATE:0:8}" - - -echo $PDATE -echo $VDATE -echo $EDATECYC -echo $EDATEDIR - -tarfile=enkfgdas.tar - - -hpssdir=/1year/NCEPDEV/emc-global/Andrew.Eichmann/HERA/scratch/${EXP}/${EDATE} - -hpsstar get ${hpssdir}/${tarfile} ./${EDATEDIR}/${EDATECYC}/gdas.t${EDATECYC}z.cnvstat.ensmean -hpsstar get ${hpssdir}/${tarfile} ./${EDATEDIR}/${EDATECYC}/gdas.t${EDATECYC}z.oznstat.ensmean -hpsstar get ${hpssdir}/${tarfile} ./${EDATEDIR}/${EDATECYC}/gdas.t${EDATECYC}z.radstat.ensmean -hpsstar get ${hpssdir}/${tarfile} ./${EDATEDIR}/${EDATECYC}/gdas.t${EDATECYC}z.abias_int.ensmean - -hpssdir=/1year/NCEPDEV/emc-global/Andrew.Eichmann/HERA/scratch/${EXP}/${PDATE} - -hpsstar get ${hpssdir}/${tarfile} ./${PDATEDIR}/${PDATECYC}/gdas.t${PDATECYC}z.atmf006.ensmean.nemsio - - - -imem=1 -while [[ $imem -le $NMEM_ENKF ]]; do - grpnum=`echo "( ($imem - 1) / ${NMEM_EARCGRP} ) + 1" | bc` - group="grp"`printf %02i $grpnum` - member="mem"`printf %03i $imem` - - echo $member - echo $group - - tarfile=enkfgdas_${group}.tar - hpssdir=/1year/NCEPDEV/emc-global/Andrew.Eichmann/HERA/scratch/${EXP}/${EDATE} - hpsstar get ${hpssdir}/${tarfile} ./${EDATEDIR}/${EDATECYC}/${member}/gdas.t${EDATECYC}z.cnvstat - hpsstar get ${hpssdir}/${tarfile} ./${EDATEDIR}/${EDATECYC}/${member}/gdas.t${EDATECYC}z.oznstat - hpsstar get ${hpssdir}/${tarfile} ./${EDATEDIR}/${EDATECYC}/${member}/gdas.t${EDATECYC}z.radstat - hpsstar get ${hpssdir}/${tarfile} ./${EDATEDIR}/${EDATECYC}/${member}/gdas.t${EDATECYC}z.atmanl.nemsio - - hpssdir=/1year/NCEPDEV/emc-global/Andrew.Eichmann/HERA/scratch/${EXP}/${PDATE} - hpsstar get ${hpssdir}/${tarfile} ./${PDATEDIR}/${PDATECYC}/${member}/gdas.t${PDATECYC}z.atmf006s.nemsio - - (( imem = $imem + 1 )) -done - - - diff --git a/util/EFSOI_Utilities/scripts/osense.py b/util/EFSOI_Utilities/scripts/osense.py deleted file mode 100644 index e4bad5028e..0000000000 --- a/util/EFSOI_Utilities/scripts/osense.py +++ /dev/null @@ -1,136 +0,0 @@ -import struct -from sys import exit, argv -import pandas as pd - -def read_osense( filename): - -#! Structure for observation sensitivity information output -#type obsense_header -# sequence -# integer(i_kind) :: idate ! Base date (initial date) -# integer(i_kind) :: obsnum ! Observation number (total) -# integer(i_kind) :: convnum ! Observation number (conventional) -# integer(i_kind) :: oznum ! Observation number (ozone) -# integer(i_kind) :: satnum ! Observation number (satellite) -# integer(i_kind) :: npred ! Number of predictors for bias correction -# integer(i_kind) :: nanals ! Number of members -#end type obsense_header -# -# where i_kind - generic specification kind for default integer -# is a long integer - - -# -#! Type definition for observation sensitivity information file -#type obsense_info -# sequence -# real(r_single) :: obfit_prior ! Observation fit to the first guess -# real(r_single) :: obsprd_prior ! Spread of observation prior -# real(r_single) :: ensmean_obnobc ! Ensemble mean first guess (no bias correction) -# real(r_single) :: ensmean_ob ! Ensemble mean first guess (bias corrected) -# real(r_single) :: ob ! Observation value -# real(r_single) :: oberrvar ! Observation error variance -# real(r_single) :: lon ! Longitude -# real(r_single) :: lat ! Latitude -# real(r_single) :: pres ! Pressure -# real(r_single) :: time ! Observation time -# real(r_single) :: oberrvar_orig ! Original error variance -# integer(i_kind) :: stattype ! Observation type -# character(len=20) :: obtype ! Observation element / Satellite name -# integer(i_kind) :: indxsat ! Satellite index (channel) set to zero -# real(r_single) :: osense_kin ! Observation sensitivity (kinetic energy) [J/kg] -# real(r_single) :: osense_dry ! Observation sensitivity (Dry total energy) [J/kg] -# real(r_single) :: osense_moist ! Observation sensitivity (Moist total energy) [J/kg] -#end type obsense_info -# - - - -#if len(sys.argv) > 0: -#if len(argv) > 0: -# file = argv[1] - - - datacolumns = [ 'obfit_prior', - 'obsprd_prior', - 'ensmean_obnobc', - 'ensmean_ob', - 'ob', - 'oberrvar', - 'lon', - 'lat', - 'pres', - 'time', - 'oberrvar_orig', - 'stattype', - 'obtype', - 'indxsat', - 'osense_kin', - 'osense_dry', - 'osense_moist' ] - - - headerf=' 0: -if len(argv) > 0: - file = argv[1] - - -datacolumns = [ 'obfit_prior', - 'obsprd_prior', - 'ensmean_obnobc', - 'ensmean_ob', - 'ob', - 'oberrvar', - 'lon', - 'lat', - 'pres', - 'time', - 'oberrvar_orig', - 'stattype', - 'obtype', - 'indxsat', - 'osense_kin', - 'osense_dry', - 'osense_moist' ] - - - -headerf=' 3176 ): -# print('indxsat=',indxsat,' ',idate) - - print('read satellite data...') - -convdata = pd.DataFrame( convdatatmp, columns = datacolumns ) - -satdata = pd.DataFrame( satdatatmp, columns = datacolumns ) - -print(convdata.head()) -print(convdata.tail()) -print(satdata.head()) -print(satdata.tail()) -#print(lon, lat, pres, time,obtype) - diff --git a/util/EFSOI_Utilities/src/CMakeLists.txt b/util/EFSOI_Utilities/src/CMakeLists.txt deleted file mode 100644 index 810ee2b36d..0000000000 --- a/util/EFSOI_Utilities/src/CMakeLists.txt +++ /dev/null @@ -1,40 +0,0 @@ -cmake_minimum_required(VERSION 2.6) -if(BUILD_EFSOI) - # enable_language (Fortran) - # set(Fortran_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include") - # set(CMAKE_Fortran_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include") - # set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") - - - file(GLOB LOCAL_SRC ${CMAKE_CURRENT_SOURCE_DIR}/*90) - -# set(Util_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include/efsoi") - - message(echo ${PROJECT_BINARY_DIR}) - - set_source_files_properties( ${LOCAL_SRC} PROPERTIES COMPILE_FLAGS ${ENKF_Fortran_FLAGS} ) - - add_executable(global_efsoi.x ${LOCAL_SRC} ) - - set_target_properties( global_efsoi.x PROPERTIES COMPILE_FLAGS ${ENKF_Fortran_FLAGS} ) - -# set_target_properties( global_efsoi.x PROPERTIES Fortran_MODULE_DIRECTORY ${Util_MODULE_DIRECTORY} ) - - SET( CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OMPFLAG}" ) - - #include_directories( "${PROJECT_BINARY_DIR}/include/global" ${CMAKE_CURRENT_BINARY_DIR} ${CORE_INCS} ${SIGIOINC} ${ENKFINC} ${FV3GFS_NCIO_INCS} ${NETCDF_INCLUDE_DIRS} ${MPI_Fortran_INCLUDE_DIRS} ${MPI_Fortran_INCLUDE_PATH} ) - -# copied from src/enkf/CMakeLists.txt - include_directories(${CMAKE_CURRENT_BINARY_DIR} "${PROJECT_BINARY_DIR}/include/wrf" "${PROJECT_BINARY_DIR}/include/global" ${CMAKE_CURRENT_BINARY_DIR}/.. ${MPI_Fortran_INCLUDE_DIRS} ${MPI_Fortran_INCLUDE_PATH} ${CORE_INCS} ${NETCDF_INCLUDE_DIRS} ${NCDIAG_INCS} ${FV3GFS_NCIO_INCS}) - -# target_link_libraries( global_efsoi.x ${GSISHAREDLIB} ${GSILIB} enkflib enkfdeplib ${NCDIAG_LIBRARIES} ${LAPACK_LIBRARIES} ${GSILIB} ${GSISHAREDLIB} ${CORE_LIBRARIES} ${CORE_BUILT} ${NETCDF_LIBRARIES_F90} ${NETCDF_LIBRARIES} ${MPI_Fortran_LIBRARIES} ${FV3GFS_NCIO_LIBRARIES} ) - -#mostly copied from src/enkf/CMakeLists.txt - target_link_libraries( global_efsoi.x enkflib enkfdeplib ${GSILIB} ${GSISHAREDLIB} ${CORE_LIBRARIES} - ${MPI_Fortran_LIBRARIES} ${LAPACK_LIBRARIES} ${NETCDF_Fortran_LIBRARIES} ${NETCDF_C_LIBRARIES} - ${FV3GFS_NCIO_LIBRARIES} ${HDF5_Fortran_HL_LIBRARIES} - ${EXTRA_LINKER_FLAGS} ${HDF5_LIBRARIES} ${CURL_LIBRARIES} ${GSI_LDFLAGS} ${CORE_BUILT} ${CORE_LIBRARIES} ${CORE_BUILT} ${NCDIAG_LIBRARIES}) -endif() - - - diff --git a/util/EFSOI_Utilities/src/CMakeLists.txtORIG b/util/EFSOI_Utilities/src/CMakeLists.txtORIG deleted file mode 100644 index 99e2c092f2..0000000000 --- a/util/EFSOI_Utilities/src/CMakeLists.txtORIG +++ /dev/null @@ -1,29 +0,0 @@ -cmake_minimum_required(VERSION 2.6) -if(BUILD_EFSOI) - # enable_language (Fortran) - # set(Fortran_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include") - # set(CMAKE_Fortran_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include") - # set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") - - - file(GLOB LOCAL_SRC ${CMAKE_CURRENT_SOURCE_DIR}/*90) - -# set(Util_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include/efsoi") - - message(echo ${PROJECT_BINARY_DIR}) - - set_source_files_properties( ${LOCAL_SRC} PROPERTIES COMPILE_FLAGS ${ENKF_Fortran_FLAGS} ) - - add_executable(global_efsoi.x ${LOCAL_SRC} ) - - set_target_properties( global_efsoi.x PROPERTIES COMPILE_FLAGS ${ENKF_Fortran_FLAGS} ) - -# set_target_properties( global_efsoi.x PROPERTIES Fortran_MODULE_DIRECTORY ${Util_MODULE_DIRECTORY} ) - - SET( CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OMPFLAG}" ) - - include_directories( "${PROJECT_BINARY_DIR}/include/global" ${CMAKE_CURRENT_BINARY_DIR}${CORE_INCS} ${SIGIOINC} ${ENKFINC} ${FV3GFS_NCIO_INCS} ) - - target_link_libraries( global_efsoi.x ${GSISHAREDLIB} ${GSILIB} enkflib enkfdeplib ${NCDIAG_LIBRARIES} ${LAPACK_LIBRARIES} ${GSILIB} ${GSISHAREDLIB} ${CORE_LIBRARIES} ${CORE_BUILT} ${NETCDF_LIBRARIES_F90} ${NETCDF_LIBRARIES} ${MPI_Fortran_LIBRARIES} ${FV3GFS_NCIO_LIBRARIES} ) - -endif() diff --git a/util/EFSOI_Utilities/src/efsoi.f90 b/util/EFSOI_Utilities/src/efsoi.f90 deleted file mode 100644 index a761b106f0..0000000000 --- a/util/EFSOI_Utilities/src/efsoi.f90 +++ /dev/null @@ -1,331 +0,0 @@ -module efsoi -!$$$ module documentation block -! -! module: efsoi Update observation impact estimates -! with the EnSRF framework. -! -! prgmmr: ota org: np23 date: 2011-12-01 -! prgmmr: groff org: emc date: 2018-05-24 -! -! abstract: -! Computes the observation impact estimates using the EnKF analysis products. -! Formulation is based on Kalnay et al (2012, Tellus A, submitted). -! Parallel processing is following the way of the serial EnSRF. -! -! Public Subroutines: -! efsoi_update: performs the Forecast Sensitivity to Observations update within -! the EnSRF framework. The EnKF other than the serial EnSRF -! (e.g. the LETKF) can be also used here (the method is -! independent on the type of EnKF). -! -! Public Variables: None -! -! Modules Used: kinds, constants, params, covlocal, mpisetup, loadbal, statevec, -! kdtree2_module, obsmod, gridinfo, obs_sensitivity -! -! program history log: -! 2011-12-01 Created from the LETKF core module. -! 2012-04-04 Changed to EnSRF framework for saving memory -! 2018-05-24 Renamed module added non-bias corrected EFSOI -! -! attributes: -! language: f95 -! -!$$$ - -use mpisetup -use covlocal, only: taper -use kinds, only: r_double, i_kind, r_kind -use kdtree2_module, only: kdtree2_r_nearest, kdtree2_result -use loadbal_efsoi, only: numptsperproc, indxproc, lnp_chunk, kdtree_grid, & - iprocob, indxob_chunk, anal_obchunk_prior, numobsperproc, & - indxproc_obs, nobs_max -! indxproc_obs, nobs_max, anal_chunk -!use statevec_efsoi, only: anal_chunk, fcerror_chunk -use scatter_chunks_efsoi, only: fcerror_chunk, anal_chunk -!use scatter_chunks_efsoi, only: fcerror_chunk -use enkf_obsmod, only: oberrvar, ob, ensmean_ob, obloc, obloclon, obloclat, oblnp, & - obtime, nobstot, corrlengthsq, lnsigl, obtimel, anal_ob -use constants, only: constants_initialized, pi, zero, one -use params, only: nanals, corrlengthnh, corrlengthtr, corrlengthsh, & - tar_minlon, tar_maxlon, tar_minlat, tar_maxlat, & - tar_minlev, tar_maxlev -use gridinfo_efsoi, only: nlevs_pres, lonsgrd, latsgrd, id_u, id_v, id_t, id_q, id_ps -use enkf_obs_sensitivity, only: obsense_kin, obsense_dry, obsense_moist, adloc_chunk - -implicit none - -private -public :: efsoi_update - -contains - -subroutine efsoi_update() -implicit none -real(r_kind),allocatable,dimension(:) :: djdy_kin, djdy_dry, djdy_moist -real(r_kind),allocatable,dimension(:) :: uvwork, tpwork, qwork -real(r_kind),allocatable,dimension(:) :: taper_disgrd -real(r_kind),dimension(nobstot) :: obdep, oberinv -real(r_single),dimension(nobstot) :: invobtimel, invlnsigl -real(r_single),dimension(nobstot) :: invcorlen -real(r_kind),allocatable,dimension(:,:) :: buffertmp -type(kdtree2_result),allocatable,dimension(:) :: sresults1 -real(r_kind) :: anal_obtmp(nanals) -real(r_kind) :: r,taper1,taper3,taperv -real(r_single) :: lnsig -real(r_double) :: t1, t2, t3, t4, tbegin, tend -integer(i_kind) :: np, nob, nob1, nob2, nobm, npob, nf2, i, ii, npt, nanal, nn -integer(i_kind) :: nnpt, nsame, nskip, ngrd1 -integer(i_kind) :: ierr -logical :: kdgrid - - -if (.not. constants_initialized) then - print *,'constants not initialized (with init_constants, init_constants_derived)' - call stop2(28) -end if - -allocate(sresults1(nlevs_pres*numptsperproc(nproc+1))) -allocate(taper_disgrd(nlevs_pres*numptsperproc(nproc+1))) - -! Compute the inverse of cut-off length and -! the observation departure from first guess -!$omp parallel do private(nob) -do nob=1,nobstot - invcorlen(nob) = one/corrlengthsq(nob) - invlnsigl(nob) = one/lnsigl(nob) - invobtimel(nob)= one/obtimel(nob) - oberinv(nob) = one/oberrvar(nob) - obdep(nob) = ob(nob)-ensmean_ob(nob) -end do - - -kdgrid=associated(kdtree_grid) -allocate(djdy_kin(nobstot)) -allocate(djdy_dry(nobstot)) -allocate(djdy_moist(nobstot)) -djdy_kin(1:nobstot) = zero -djdy_dry(1:nobstot) = zero -djdy_moist(1:nobstot) = zero -allocate(uvwork(nanals)) -allocate(tpwork(nanals)) -allocate(qwork(nanals)) - -nsame = 0 -nskip = 0 -nobm = 1 -t2 = zero -t3 = zero -t4 = zero -tbegin = mpi_wtime() - -! -! Loop for each observations -obsloop: do nob=1,nobstot - - - if (nproc == 0 .and. modulo(nob , 10000) == 0 ) print *,'doing nob=',nob - - t1 = mpi_wtime() - - if(oberrvar(nob) > 1.e10_r_kind .or. abs(obtime(nob)) >= obtimel(nob))then - nskip = nskip + 1 - cycle obsloop - end if - - ! what processor is this ob on? - npob = iprocob(nob) - - ! get anal_ob from that processor and send to other processors. - if (nproc == npob) then - nob1 = indxob_chunk(nob) - anal_obtmp(1:nanals) = anal_obchunk_prior(:,nob1) - end if - - call mpi_bcast(anal_obtmp,nanals,mpi_realkind,npob,mpi_comm_world,ierr) - - anal_obtmp(1:nanals) = anal_obtmp(1:nanals) * oberinv(nob) - - t2 = t2 + mpi_wtime() - t1 - t1 = mpi_wtime() - - - - - - ! ============================================================ - ! Only need to recalculate nearest points when lat/lon is different - ! ============================================================ - if(nob == 1 .or. & - abs(obloclat(nob)-obloclat(nobm)) .gt. tiny(obloclat(nob)) .or. & - abs(obloclon(nob)-obloclon(nobm)) .gt. tiny(obloclon(nob)) .or. & - abs(corrlengthsq(nob)-corrlengthsq(nobm)) .gt. tiny(corrlengthsq(nob))) then - - nobm=nob - ! determine localization length scales based on latitude of ob. - - nf2=0 - ! search analysis grid points for those within corrlength of - ! ob being assimilated (using a kd-tree for speed). - if (kdgrid) then - call kdtree2_r_nearest(tp=kdtree_grid,qv=obloc(:,nob), & - r2=corrlengthsq(nob), & - nfound=nf2,nalloc=nlevs_pres*numptsperproc(nproc+1), & - results=sresults1) - else - ! use brute force search if number of grid points on this proc <= 3 - do nn=1,nlevs_pres - do npt=1,numptsperproc(nproc+1) - nnpt = nlevs_pres * (npt-1) + nn - r = sum( (obloc(:,nob)-adloc_chunk(:,nnpt))**2, 1 ) - if (r < corrlengthsq(nob)) then - nf2 = nf2 + 1 - sresults1(nf2)%idx = nnpt - sresults1(nf2)%dis = r - end if - end do - end do - end if - - !$omp parallel do private(nob1) - do nob1=1,nf2 - taper_disgrd(nob1) = taper(sqrt(sresults1(nob1)%dis*invcorlen(nob))) - end do - - else - nsame=nsame+1 - end if - - t3 = t3 + mpi_wtime() - t1 - t1 = mpi_wtime() - - ! ============================================================ - ! forecast error sensitivity to the observations - ! ============================================================ - if (nf2 > 0) then - - taper3=taper(obtime(nob)*invobtimel(nob)) - - uvwork(1:nanals) = zero - tpwork(1:nanals) = zero - qwork(1:nanals) = zero - - ! rho * X_f^T (e_t|0 + e_t|-6) / 2 - ! contributions from (U,V), (T,Ps) and Q are computed separately - - do ii=1,nf2 - taper1=taper_disgrd(ii)*taper3 - - i = (sresults1(ii)%idx-1)/nlevs_pres+1 - - ngrd1=indxproc(nproc+1,i) - - if(tar_minlon <= tar_maxlon .and. & - & (lonsgrd(ngrd1) < tar_minlon .or. lonsgrd(ngrd1) > tar_maxlon .or. & - & latsgrd(ngrd1) < tar_minlat .or. latsgrd(ngrd1) > tar_maxlat)) & - & cycle - - if(tar_minlon > tar_maxlon .and. & - & ((lonsgrd(ngrd1) < tar_minlon .and. lonsgrd(ngrd1) > tar_maxlon) .or. & - & latsgrd(ngrd1) < tar_minlat .or. latsgrd(ngrd1) > tar_maxlat)) & - & cycle - - nn = sresults1(ii)%idx - (i-1)*nlevs_pres - - if((tar_minlev /= 1 .or. nn /= nlevs_pres) & - & .and. (nn > tar_maxlev .or. nn < tar_minlev)) cycle - lnsig = abs(lnp_chunk(i,nn)-oblnp(nob)) - - if(lnsig < lnsigl(nob))then - taperv=taper1*taper(lnsig*invlnsigl(nob)) - else - cycle - end if - - if(nn == nlevs_pres) then - do nanal=1,nanals - tpwork(nanal) = tpwork(nanal) + taperv & - & * anal_chunk(nanal,i,id_ps) * fcerror_chunk(i,id_ps) - end do - cycle - end if - - - ! - do nanal=1,nanals - uvwork(nanal) = uvwork(nanal) + taperv & - & * (anal_chunk(nanal,i,id_u(nn)) * fcerror_chunk(i,id_u(nn)) & - & + anal_chunk(nanal,i,id_v(nn)) * fcerror_chunk(i,id_v(nn))) - - tpwork(nanal) = tpwork(nanal) + taperv & - & * anal_chunk(nanal,i,id_t(nn)) * fcerror_chunk(i,id_t(nn)) - - if(id_q(nn) > 0) qwork(nanal) = qwork(nanal) + taperv & - & * anal_chunk(nanal,i,id_q(nn)) * fcerror_chunk(i,id_q(nn)) - end do - - end do - - ! R^-1 HX_a [rho * X_f^T (e_t|0 + e_t|-6) / 2] - do nanal=1,nanals - djdy_kin(nob) = djdy_kin(nob) + anal_obtmp(nanal) * uvwork(nanal) - djdy_dry(nob) = djdy_dry(nob) + anal_obtmp(nanal) * tpwork(nanal) - djdy_moist(nob) = djdy_moist(nob) + anal_obtmp(nanal) * qwork(nanal) - end do - - end if ! end of, if (nf2 > 0) then - - - t4 = t4 + mpi_wtime() - t1 - t1 = mpi_wtime() - -end do obsloop - -tend = mpi_wtime() -if (nproc == 0 .or. nproc == numproc-1) print *,'time to process FSO on gridpoint = ',tend-tbegin,t2,t3,t4,' secs' - -if (nproc == 0 .and. nskip > 0) print *,nskip,' out of',nobstot,'obs skipped' -if (nproc == 0 .and. nsame > 0) print *,nsame,' out of', nobstot-nskip,' same lat/lon' - -! Observation sensitivity post process -!$omp parallel do private(nob) -do nob=1,nobstot - djdy_dry(nob) = djdy_dry(nob) + djdy_kin(nob) - djdy_moist(nob) = djdy_moist(nob) + djdy_dry(nob) - obsense_kin(nob) = djdy_kin(nob) * obdep(nob) - obsense_dry(nob) = djdy_dry(nob) * obdep(nob) - obsense_moist(nob) = djdy_moist(nob) * obdep(nob) -end do - -! Gathering analysis perturbations projected on the observation space -if(nproc /= 0) then - call mpi_send(anal_obchunk_prior,numobsperproc(nproc+1)*nanals,mpi_real4,0, & - 1,mpi_comm_world,ierr) -else - allocate(anal_ob(1:nanals,nobstot)) - allocate(buffertmp(nanals,nobs_max)) - do np=1,numproc-1 - call mpi_recv(buffertmp,numobsperproc(np+1)*nanals,mpi_real4,np, & - 1,mpi_comm_world,mpi_status,ierr) - do nob1=1,numobsperproc(np+1) - nob2 = indxproc_obs(np+1,nob1) - anal_ob(:,nob2) = buffertmp(:,nob1) - end do - end do - do nob1=1,numobsperproc(1) - nob2 = indxproc_obs(1,nob1) - anal_ob(:,nob2) = anal_obchunk_prior(:,nob1) - end do - deallocate(buffertmp) -end if - -t1 = mpi_wtime() - tend -if (nproc == 0) print *,'time to observation sensitivity post process = ',t1,' secs' - -deallocate(anal_obchunk_prior) -deallocate(sresults1) -deallocate(djdy_kin,djdy_dry,djdy_moist,taper_disgrd,uvwork,tpwork,qwork) - -return -end subroutine efsoi_update -end module efsoi diff --git a/util/EFSOI_Utilities/src/efsoi_main.f90 b/util/EFSOI_Utilities/src/efsoi_main.f90 deleted file mode 100644 index 008f6b4231..0000000000 --- a/util/EFSOI_Utilities/src/efsoi_main.f90 +++ /dev/null @@ -1,166 +0,0 @@ -program efsoi_main -!$$$ main program documentation block -! -! program: efsoi_main high level driver program for -! efsoi calculations. -! -! prgmmr: Ota org: EMC/JMA date: 2012 -! prgmmr: Groff org: EMC date: 2018 -! -! abstract: This is the main program for EFSOI code. It does the following: -! a) initialize MPI, read EFSOI namelist from efsoi.nml on each task. -! b) reads observation sensitivity files. -! c) read horizontal grid information (lat/lon of each grid point) and -! pressure at each grid point/vertical level. -! d) decomposition of horizontal grid points and observation -! priors to minimize load imbalance for EFSOI calcs. -! e) read forecast and analysis states necessary for -! EFSOI calculations. -! f) Initialize/allocate for EFSOI moist/dry/kinetic -! total output -! g) Estimate the location of observation response -! (i.e. advect the localization). Default estimation -! approach is to multiply meridional and zonal wind -! by 0.75. -! h) Perform EFSOI calculations for all observations -! considered for assimilation during the EnSRF update. -! i) write EFSOI calculations and ancillary EFSOI -! information to file -! j) deallocate all allocatable arrays, finalize MPI. -! -! program history log: -! 2018-04-30 Initial development. Adaptation of enkf_main program -! towards EFSOI calculation process -! -! usage: -! input files: Update to FV3 nomenclature -! sigfAT_YYYYMMDDHH_mem* - Ensemble of forecasts valid at advance time -! (AT) hours beyond the initial time -! sigfAT_YYYYMMDDHH_ensmean - Ensemble mean forecast AT hours from -! initial time -! sigfAT+6_YYYYMMDDHH-6_ensmean - Ensemble mean forecast AT+6 hours -! from initial time minus 6 hours -! siganl.YYYYMMDDHH.gdas - -! output files: -! obimpact_YYYYMMDDHH.dat - observation impact file -! -! comments: This program is a wrapper for components needed to peform -! EFSOI calculations -! -! attributes: -! language: f95 -! -!$$$ - - use kinds, only: r_double,i_kind - ! reads namelist parameters. - ! applying enkf namelist apparatus - use params, only : read_namelist,nanals - ! mpi functions and variables. - use mpisetup, only: mpi_initialize, mpi_initialize_io, mpi_cleanup, nproc, & - mpi_wtime, mpi_comm_world - ! model state vector - use statevec_efsoi, only: read_state_efsoi, statevec_cleanup_efsoi, init_statevec_efsoi - ! load balancing - use loadbal_efsoi, only: load_balance_efsoi, loadbal_cleanup_efsoi - ! efsoi update - use efsoi, only: efsoi_update - ! Observation sensitivity usage - use enkf_obs_sensitivity, only: init_ob_sens, print_ob_sens, destroy_ob_sens, read_ob_sens - use loc_advection, only: loc_advection_efsoi - ! Scatter chunks for EFSOI - use scatter_chunks_efsoi, only: scatter_chunks_ob_impact - use omp_lib, only: omp_get_max_threads - ! Needed for rad fixed files - use radinfo, only: init_rad, init_rad_vars - - implicit none - integer :: ierr, nth - real(r_double) t1,t2 - - ! 0. initialize MPI. - call mpi_initialize() - if (nproc==0) call w3tagb('EFSOI_CALC',2018,0319,0055,'NP25') - - call init_rad() - - ! 1. read namelist. - call read_namelist() - - ! 2. initialize MPI communicator for IO tasks. - call mpi_initialize_io(nanals) - - ! 3. Initialize rad vars - call init_rad_vars() - - ! 4. read the necessary inputs for the EFSOI calculation from file - t1 = mpi_wtime() - call read_ob_sens() - t2 = mpi_wtime() - if (nproc == 0) print *, 'time in read_ob_sens = ',t2-t1,'on proc', nproc - - nth= omp_get_max_threads() - if(nproc== 0)write(6,*) 'efsoi_main: number of threads ',nth - - ! 5. Halt processors until all are completed - call mpi_barrier(mpi_comm_world, ierr) - - ! 6. Initialize state vector information - call init_statevec_efsoi() - - ! 7. read in ensemble forecast members, - ! valid at the evaluation forecast - ! time, distribute pieces to each task. - t1 = mpi_wtime() - call read_state_efsoi() - t2 = mpi_wtime() - if (nproc == 0) print *,'time in read_state_efsoi =',t2-t1,'on proc',nproc - - ! 8. do load balancing (partitioning of grid points - ! and observations among processors) - t1 = mpi_wtime() - call load_balance_efsoi() - t2 = mpi_wtime() - if (nproc == 0) print *,'time in load_balance_efsoi =',t2-t1,'on proc',nproc - - ! 9. apply scattering of efsoi chunks - t1 = mpi_wtime() - call scatter_chunks_ob_impact() ! ensemble scattering - t2 = mpi_wtime() - if (nproc == 0) print *,'time to scatter observation impact chunks =',t2-t1,'on proc',nproc - - ! 10. Initialize EFSOI variables - t1 = mpi_wtime() - call init_ob_sens() - t2 = mpi_wtime() - if (nproc == 0) print *,'time to allocate ob sensitivity variables =',t2-t1,'on proc',nproc - - ! 11. Calculate the estimated location of observation - ! response at evaluation time - t1 = mpi_wtime() - call loc_advection_efsoi() - t2 = mpi_wtime() - if (nproc == 0) print *,'time in loc_advection_efsoi =',t2-t1,'on proc',nproc - - ! 12. Perform the EFSOI calcs - t1 = mpi_wtime() - call efsoi_update() - t2 = mpi_wtime() - if (nproc == 0) print *,'time in efsoi_update =',t2-t1,'on proc',nproc - - ! 13. print EFSOI sensitivity i/o on root task. - t1 = mpi_wtime() - call print_ob_sens() - t2 = mpi_wtime() - if (nproc == 0) print *,'time needed to write observation impact file =',t2-t1,'on proc',nproc - - ! 14. Cleanup for EFSOI configuration - call statevec_cleanup_efsoi() - call loadbal_cleanup_efsoi() - call destroy_ob_sens() - - ! 15. finalize MPI. - if (nproc==0) call w3tage('EFSOI_CALC') - call mpi_cleanup() - -end program efsoi_main diff --git a/util/EFSOI_Utilities/src/gridinfo_efsoi.f90 b/util/EFSOI_Utilities/src/gridinfo_efsoi.f90 deleted file mode 100644 index 74d2748271..0000000000 --- a/util/EFSOI_Utilities/src/gridinfo_efsoi.f90 +++ /dev/null @@ -1,461 +0,0 @@ -module gridinfo_efsoi -!$$$ module documentation block -! -! module: gridinfo_efsoi read horizontal (lons, lats) and -! vertical (pressure) information from -! ensemble mean first guess file. -! -! prgmmr: whitaker org: esrl/psd date: 2009-02-23 -! prgmmr: groff org: emc date: 2018-05-24 -! -! abstract: This module reads gfg_YYYYMMDDHH_fhr06_ensmean, and -! extracts information about the analysis grid, including the -! longitudes and latitudes of the analysis grid points and -! the pressure on each grid point/vertical level. -! -! Public Subroutines: -! getgridinfo_efsoi: read latitudes, longitudes, pressures and orography for analysis grid, -! broadcast to each task. Compute spherical cartesian coordinate values -! for each analysis horizontal grid point. -! gridinfo_cleanup_efsoi: deallocate allocated module variables. -! -! Public Variables: -! npts: number of analysis grid points in the horizontal (from module params). -! nlevs: number of analysis vertical levels (from module params). -! ntrac: number of 'tracer' model state variables (3 for GFS, -! specific humidity, ozone and cloud condensate). -! ptop: (real scalar) pressure (hPa) at top model layer interface. -! lonsgrd(npts): real array of analysis grid longitudes (radians). -! latsgrd(npts): real array of analysis grid latitudes (radians). -! logp(npts,ndim): -log(press) for all 2d analysis grids. Assumed invariant -! in assimilation window, computed fro ensemble mean at middle of window. -! gridloc(3,npts): spherical cartesian coordinates (x,y,z) for analysis grid. -! -! Modules Used: mpisetup, params, kinds -! -! program history log: -! 2009-02-23 Initial version. -! 2016-05-02: shlyaeva: Modification for reading state vector from table -! 2016-04-20 Modify to handle the updated nemsio sig file (P, DP & DPDT removed) -! 2018-05-24 Groff: Adapting EnKF gridinfo_efsoi module to EFSOI needs. -! -! attributes: -! language: f95 -! -!$$$ - -! ========================================= -!use mpisetup, only: nproc, mpi_integer, mpi_real4, mpi_comm_world -!use params, only: datapath,nlevs,nlons,nlats,fgfileprefixes,nvars - -use mpisetup, only: nproc, mpi_integer, mpi_real4 -use mpimod, only: mpi_comm_world -use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio,use_gfs_ncio,fgfileprefixes -! ========================================= - -use kinds, only: r_kind, i_kind, r_double, r_single -use constants, only: one,zero,pi,cp,rd,grav,rearth,max_varname_length -use specmod, only: sptezv_s, sptez_s, init_spec_vars, isinitialized, asin_gaulats, & - ndimspec => nc -use reducedgrid_mod, only: reducedgrid_init, regtoreduced, reducedtoreg,& - nptsred, lonsred, latsred -use mpeu_util, only: getindex - -!use statevec_efsoi, only: clevels !, cvars2d, cvars3d, nc3d -implicit none -private -public :: getgridinfo_efsoi, gridinfo_cleanup_efsoi -integer(i_kind),public :: nlevs_pres, idvc -real(r_single),public :: ptop -real(r_single),public, allocatable, dimension(:) :: lonsgrd, latsgrd -! arrays passed to kdtree2 routines must be single -real(r_single),public, allocatable, dimension(:,:) :: gridloc -real(r_single),public, allocatable, dimension(:,:) :: logp -integer,public :: npts -integer,public :: ntrunc -integer(i_kind),public, allocatable, dimension(:) :: id_u, id_v, id_t, id_q -integer(i_kind),public :: id_ps, ncdim -! supported variable names in anavinfo -character(len=max_varname_length),public, dimension(10) :: vars3d_supported = (/'u ', 'v ', 'tv ', 'q ', 'oz ', 'cw ', 'tsen', 'prse', 'ql ', 'qi '/) -character(len=max_varname_length),public, dimension(3) :: vars2d_supported = (/'ps ', 'pst', 'sst' /) -! supported variable names in anavinfo -contains - - - - - -subroutine getgridinfo_efsoi(fileprefix, reducedgrid, clevels, nc3d) -! read latitudes, longitudes and pressures for analysis grid, -! broadcast to each task. -use sigio_module, only: sigio_head, sigio_data, sigio_sclose, sigio_sropen, & - sigio_srohdc, sigio_sclose, sigio_srhead, sigio_axdata - -use nemsio_module, only: nemsio_gfile,nemsio_open,nemsio_close,& - nemsio_getfilehead,nemsio_getheadvar,& - nemsio_readrecv,nemsio_init, nemsio_realkind - -use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& - read_attribute, close_dataset, get_dim, read_vardata - - -implicit none - -type(Dataset) :: dset -type(Dimension) :: londim,latdim,levdim -character(len=120), intent(in) :: fileprefix -logical, intent(in) :: reducedgrid - -integer(i_kind) nlevsin, ierr, iunit, k, nn, idvc -character(len=500) :: filename -real(r_kind), allocatable, dimension(:) :: ak,bk,spressmn,tmpspec -real(r_kind), allocatable, dimension(:,:) :: pressimn,presslmn,values_2d -real(r_single),allocatable,dimension(:,:,:) :: nems_vcoord -real(r_kind) kap,kapr,kap1 -real(nemsio_realkind), dimension(nlons*nlats) :: nems_wrk -type(sigio_data) sigdata -type(sigio_head) sighead -type(nemsio_gfile) :: gfile - -integer(i_kind) iret,i,j,nlonsin,nlatsin,u_ind,v_ind,tv_ind,q_ind,ps_ind -integer, intent(in) :: nc3d -integer, dimension(0:nc3d), intent(in) :: clevels - - -iunit = 77 -kap = rd/cp -kapr = cp/rd -kap1 = kap + one -nlevs_pres=nlevs+1 - - -if (nproc .eq. 0) then - filename = trim(adjustl(datapath))//trim(adjustl(fileprefix))//"ensmean" - print *,'reading file: ',trim(filename) - - if (use_gfs_nemsio) then - call nemsio_init(iret=iret) - if(iret/=0) then - write(6,*)'grdinfo: gfs model: problem with nemsio_init, iret=',iret, ', file: ', trim(filename) - call stop2(23) - end if - call nemsio_open(gfile,filename,'READ',iret=iret) - if (iret/=0) then - write(6,*)'grdinfo: gfs model: problem with nemsio_open, iret=',iret, ', file: ', trim(filename) - call stop2(23) - endif - call nemsio_getfilehead(gfile,iret=iret, dimx=nlonsin, dimy=nlatsin,& - dimz=nlevsin,jcap=ntrunc,idvc=idvc) - ! set ntrunc to nlats if missing - ! (only used for inflation smoothing and mass balance adjustment if use_gfsnemsio = T) - ! FV3GFS write component does not include JCAP, infer from nlatsin - if (ntrunc < 0) ntrunc = nlatsin-2 - if (iret/=0) then - write(6,*)'grdinfo: gfs model: problem with nemsio_getfilehead, iret=',iret, ', file: ', trim(filename) - call stop2(23) - endif - print *,'ntrunc = ',ntrunc - if (nlons /= nlonsin .or. nlats /= nlatsin .or. nlevs /= nlevsin) then - print *,'incorrect dims in nemsio file' - print *,'expected',nlons,nlats,nlevs - print *,'got',nlonsin,nlatsin,nlevsin - call stop2(23) - end if - else if (use_gfs_ncio) then - - dset = open_dataset(filename) - - londim = get_dim(dset,'grid_xt'); nlonsin = londim%len - latdim = get_dim(dset,'grid_yt'); nlatsin = latdim%len - levdim = get_dim(dset,'pfull'); nlevsin = levdim%len - idvc = 2; ntrunc = nlatsin-2 - - if (nlons /= nlonsin .or. nlats /= nlatsin .or. nlevs /= nlevsin) then - print *,'incorrect dims in netcdf file' - print *,'expected',nlons,nlats,nlevs - print *,'got',nlonsin,nlatsin,nlevsin - call stop2(23) - end if - - else - ! define sighead on all tasks. - call sigio_sropen(iunit,trim(filename),iret) - if (iret /= 0) then - print *,'error reading file in gridinfo ',trim(filename),' on task',nproc - call stop2(24) - end if - call sigio_srhead(iunit,sighead,iret) - if (iret /= 0) then - print *,'error reading file in gridinfo ',trim(filename),' on task',nproc - call stop2(24) - end if - call sigio_sclose(iunit,iret) - ntrunc = sighead%jcap - endif -endif -! ======================================== - - -call mpi_bcast(ntrunc,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) - -! initialize spectral module on all tasks. -if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4) - -if (nproc .eq. 0) then - ! get pressure, lat/lon information from ensemble mean file. - allocate(presslmn(nlons*nlats,nlevs)) - allocate(pressimn(nlons*nlats,nlevs+1)) - allocate(spressmn(nlons*nlats)) - - ! --- nemsio data ------------------------------ - if (use_gfs_nemsio) then - call nemsio_readrecv(gfile,'pres','sfc',1,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'grdinfo: gfs model: problem with nemsio_readrecv(ps), iret=',iret - call stop2(23) - endif - - ! Extract vertical coordinate descriptions nems_vcoord. - ! nems_vcoord(gfshead%levs+1,3,2) dimension is hardwired here. - ! Present NEMSIO modules do not allow flexibility of 2nd and 3rd - ! array dimension for nems_vcoord, for now, it is hardwired as - ! (levs,3,2) If NEMS changes the setting of vcoord dimension, - ! GSI needs to update its setting of nems_vcoord accordingly. - - if (allocated(nems_vcoord)) deallocate(nems_vcoord) - allocate(nems_vcoord(nlevs_pres,3,2)) - call nemsio_getfilehead(gfile,iret=iret,vcoord=nems_vcoord) - if ( iret /= 0 ) then - write(6,*)' gridinfo_efsoi: ***ERROR*** problem reading header ', & - 'vcoord, Status = ',iret - call stop2(99) - endif - - spressmn = 0.01_r_kind*nems_wrk ! convert ps to millibars. - !print *,'min/max spressmn = ',minval(spressmn),maxval(spressmn) - - allocate(ak(nlevs+1),bk(nlevs+1)) - - if ( idvc == 0 ) then ! sigma coordinate, old file format. - ak = zero - bk = nems_vcoord(1:nlevs+1,1,1) - elseif ( idvc == 1 ) then ! sigma coordinate - ak = zero - bk = nems_vcoord(1:nlevs+1,2,1) - elseif ( idvc == 2 .or. idvc == 3 ) then ! hybrid coordinate - ak = 0.01_r_kind*nems_vcoord(1:nlevs+1,1,1) ! convert to mb - bk = nems_vcoord(1:nlevs+1,2,1) - else - write(6,*)'gridinfo_efsoi: ***ERROR*** INVALID value for idvc=',idvc - call stop2(85) - endif - - ! pressure at interfaces - do k=1,nlevs+1 - pressimn(:,k) = ak(k)+bk(k)*spressmn(:) - enddo - call nemsio_close(gfile, iret=iret) - ptop = ak(nlevs+1) - deallocate(ak,bk) - - ! --- NetCDF data ------------------------------ - else if (use_gfs_ncio) then - call read_vardata(dset, 'pressfc', values_2d,errcode=iret) - if (iret /= 0) then - print *,'error reading ps in gridinfo_gfs' - call stop2(11) - endif - ! convert to 1d array, units to millibars, flip so lats go N to S. - spressmn = 0.01_r_kind*reshape(values_2d,(/nlons*nlats/)) - call read_attribute(dset, 'ak', ak) - call read_attribute(dset, 'bk', bk) - - call close_dataset(dset) - - ! pressure at interfaces - do k=1,nlevs+1 - pressimn(:,k) = 0.01_r_kind*ak(nlevs-k+2)+bk(nlevs-k+2)*spressmn(:) - enddo - ptop = 0.01_r_kind*ak(1) - deallocate(ak,bk,values_2d) - - ! --- other data ------------------------------ - else - - ! get pressure from ensemble mean, - ! distribute to all processors. - call sigio_srohdc(iunit,trim(filename), & - sighead,sigdata,iret) - if (iret /= 0) then - print *,'error reading file in gridinfo',trim(filename) - call stop2(24) - end if - nlevsin = sighead%levs - if (nlevs .ne. nlevsin) then - print *,'error reading input file in gridinfo - nlevs != ',nlevsin,nlevs - call stop2(24) - end if - allocate(ak(nlevs+1),bk(nlevs+1)) - if (sighead%idvc == 0) then ! sigma coordinate, old file format. - ak = zero - bk = sighead%si(1:nlevs+1) - else if (sighead%idvc == 1) then ! sigma coordinate - ak = zero - bk = sighead%vcoord(1:nlevs+1,2) - else if (sighead%idvc == 2 .or. sighead%idvc == 3) then ! hybrid coordinate - ak = 0.01_r_kind*sighead%vcoord(1:nlevs+1,1) ! convert to mb - bk = sighead%vcoord(1:nlevs+1,2) - else - print *,'unknown vertical coordinate type',sighead%idvc - call stop2(24) - end if - allocate(tmpspec(ndimspec)) - tmpspec = sigdata%ps - call sptez_s(tmpspec,spressmn,1) - deallocate(tmpspec) - spressmn = 10._r_kind*exp(spressmn) - ! pressure at interfaces - do k=1,nlevs+1 - pressimn(:,k) = ak(k)+bk(k)*spressmn(:) - enddo - call sigio_axdata(sigdata,iret) - ptop = ak(nlevs+1) - deallocate(ak,bk) - endif - - if (reducedgrid) then - call reducedgrid_init(nlons,nlats,asin_gaulats) - npts = nptsred - else - npts = nlons*nlats - end if - - allocate(latsgrd(npts),lonsgrd(npts)) - allocate(logp(npts,nlevs_pres)) ! log(ens mean first guess press) on mid-layers - allocate(gridloc(3,npts)) - - - !==> pressure at interfaces. - if (reducedgrid) then - lonsgrd(:) = lonsred(:) - latsgrd(:) = latsred(:) - else - nn = 0 - do j=1,nlats - do i=1,nlons - nn = nn + 1 - lonsgrd(nn) = 2._r_single*pi*float(i-1)/nlons - latsgrd(nn) = asin_gaulats(j) - enddo - enddo - endif - - - - do k=1,nlevs - ! layer pressure from Phillips vertical interpolation. - presslmn(:,k) = ((pressimn(:,k)**kap1-pressimn(:,k+1)**kap1)/& - (kap1*(pressimn(:,k)-pressimn(:,k+1))))**kapr - end do - - - - print *,'ensemble mean first guess surface pressure:' - print *,minval(spressmn),maxval(spressmn) - !do k=1,nlevs - ! print *,'min/max ens mean press level',& - ! k,'=',minval(presslmn(:,k)),maxval(presslmn(:,k)) - ! print *,'min/max ens mean press interface',& - ! k,'=',minval(pressimn(:,k)),maxval(pressimn(:,k)) - !enddo - ! logp holds log(pressure) or pseudo-height on grid, for each level/variable. - - - do k=1,nlevs - ! all variables to be updated are on mid-layers, not layer interfaces. - if (reducedgrid) then - call regtoreduced(presslmn(:,k),logp(:,k)) - logp(:,k) = -log(logp(:,k)) - else - logp(:,k) = -log(presslmn(:,k)) - endif - !print *,'min/max presslmn',k,minval(presslmn(:,k)),maxval(presslmn(:,k)),minval(logp(:,k)),maxval(logp(:,k)) - end do - - - if (reducedgrid) then - call regtoreduced(spressmn,logp(:,nlevs_pres)) - logp(:,nlevs_pres) = -log(logp(:,nlevs_pres)) - else - logp(:,nlevs_pres) = -log(spressmn(:)) - endif - deallocate(spressmn,presslmn,pressimn) - -end if - -call mpi_bcast(npts,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) -if (nproc .ne. 0) then - ! allocate arrays on other (non-root) tasks - allocate(latsgrd(npts),lonsgrd(npts)) - allocate(logp(npts,nlevs_pres)) ! log(ens mean first guess press) on mid-layers - allocate(gridloc(3,npts)) - ! initialize reducedgrid_mod on other tasks. - if (reducedgrid) then - call reducedgrid_init(nlons,nlats,asin_gaulats) - end if -endif - -!call mpi_bcast(logp,npts*nlevs_pres,mpi_real4,0,MPI_COMM_WORLD,ierr) -do k=1,nlevs_pres - call mpi_bcast(logp(1,k),npts,mpi_real4,0,MPI_COMM_WORLD,ierr) -enddo - -call mpi_bcast(lonsgrd,npts,mpi_real4,0,MPI_COMM_WORLD,ierr) -call mpi_bcast(latsgrd,npts,mpi_real4,0,MPI_COMM_WORLD,ierr) -call mpi_bcast(ptop,1,mpi_real4,0,MPI_COMM_WORLD,ierr) - -!==> precompute cartesian coords of analysis grid points. -do nn=1,npts - gridloc(1,nn) = cos(latsgrd(nn))*cos(lonsgrd(nn)) - gridloc(2,nn) = cos(latsgrd(nn))*sin(lonsgrd(nn)) - gridloc(3,nn) = sin(latsgrd(nn)) -end do - -! Identify EFSOI relevant state variable indices -u_ind = getindex(vars3d_supported, 'u') !< indices in the state var arrays -v_ind = getindex(vars3d_supported, 'v') ! U and V (3D) -tv_ind = getindex(vars3d_supported, 'tv') ! Tv (3D) -q_ind = getindex(vars3d_supported, 'q') ! Q (3D) -ps_ind = getindex(vars2d_supported, 'ps') ! Ps (2D) - -! Index of each elements -allocate(id_u(nlevs)) -allocate(id_v(nlevs)) -allocate(id_t(nlevs)) -allocate(id_q(nlevs)) -do k=1,nlevs - id_u(k) = clevels(u_ind-1) + k - id_v(k) = clevels(v_ind-1) + k - id_t(k) = clevels(tv_ind-1) + k - id_q(k) = clevels(q_ind-1) + k -end do -id_ps = clevels(nc3d) + ps_ind - - - - - -end subroutine getgridinfo_efsoi - -subroutine gridinfo_cleanup_efsoi() -if (allocated(lonsgrd)) deallocate(lonsgrd) -if (allocated(latsgrd)) deallocate(latsgrd) -if (allocated(logp)) deallocate(logp) -if (allocated(gridloc)) deallocate(gridloc) -if (allocated(id_u)) deallocate(id_u) -if (allocated(id_v)) deallocate(id_v) -if (allocated(id_t)) deallocate(id_t) -if (allocated(id_q)) deallocate(id_q) -end subroutine gridinfo_cleanup_efsoi - -end module gridinfo_efsoi diff --git a/util/EFSOI_Utilities/src/gridio_efsoi.f90 b/util/EFSOI_Utilities/src/gridio_efsoi.f90 deleted file mode 100644 index cd2263f6d3..0000000000 --- a/util/EFSOI_Utilities/src/gridio_efsoi.f90 +++ /dev/null @@ -1,927 +0,0 @@ - module gridio_efsoi -!$$$ module documentation block -! -! module: gridio_efsoi subroutines for reading and writing -! ensemble members files using -! EnKF internal format. A separate -! program must be run before and -! after the EnKF analysis to convert -! to and from the native model format. -! -! prgmmr: whitaker org: esrl/psd date: 2009-02-23 -! prgmmr: groff org: emc date: 2018-05-24 -! -! abstract: I/O for ensemble member files. -! -! Public Functions: -! readgriddata_efsoi -! -! this version reads nemsio files for EFSOI applications. -! -! Public Variables: None -! -! Modules Used: constants (must be pre-initialized). -! -! program history log: -! 2009-02-23 Initial version. -! 2015-06-29 Add ability to read/write multiple time levels -! 2016-05-02 shlyaeva: Modification for reading state vector from table -! 2016-04-20 Modify to handle the updated nemsio sig file (P, DP, DPDT -! removed) -! 2016-11-29 shlyaeva: Add reading/calculating tsen, qi, ql. Pass filenames and -! hours to read routine to read separately state and control files. -! Pass levels and dimenstions to read/write routines (dealing with -! prse: nlevs + 1 levels). Pass "reducedgrid" parameter. -! 2017-06-14 Adding functionality to optionally write non-inflated ensembles, -! a required input for EFSO calculations -! 2018-05-24 Pruning available EnKF nemsio read functionality for EFSOI -! application. Add additional routines to compute EFSOI relevant -! quantities from files read. -! 2019-03-13 Add precipitation components -! 2019-07-10 Add convective clouds -! -! attributes: -! language: f95 -! -!$$$ - use constants, only: zero,one,cp,fv,rd,tiny_r_kind,max_varname_length,t0c,r0_05,constants_initialized, & - tref,pref,hvap - - use params, only: nlons,nlats,nlevs,use_gfs_nemsio,pseudo_rh, & - cliptracers,datapath,imp_physics,use_gfs_ncio,cnvw_option, & - nanals, & - wmoist,andataname, & - nvars,forecast_impact,read_member_forecasts, & - read_ensmean_forecast, read_analysis_mean, & - read_member_forecasts, read_verification, & - read_member_analyses - use kinds, only: i_kind,r_double,r_kind,r_single - use gridinfo_efsoi, only: ntrunc,npts,ncdim ! getgridinfo_efsoi must be called first! - - use specmod, only: sptezv_s, sptez_s, init_spec_vars, ndimspec => nc, & - isinitialized, & - gaulats - - use reducedgrid_mod, only: regtoreduced, reducedtoreg, & - lonsperlat, nlonsfull - - ! === updated by LL on 2020.07.15 ====== - use mpisetup, only: nproc, numproc - !use mpisetup, only: nproc, mpi_integer, mpi_real4 - - use mpimod, only: mpi_comm_world, mpi_sum, mpi_real4, mpi_real8, mpi_rtype - !use mpimod, only: mpi_comm_world - ! ============================================ - - use mpeu_util, only: getindex - use nemsio_module - use loadbal_efsoi, only: numptsperproc, indxproc, npts_max - - implicit none - private - public :: readgriddata_efsoi - public :: get_weight, destroy_weight, divide_weight - real(r_kind), allocatable, dimension(:,:), save :: weight - real(r_kind), allocatable, dimension(:), save :: grweight - contains - - ! ===================================================================== - ! ===================================================================== - ! ===================================================================== - subroutine get_weight() - - use sigio_module, only: sigio_head, sigio_data, sigio_sclose, sigio_sropen, & - sigio_srohdc, sigio_sclose, sigio_aldata, sigio_axdata - use nemsio_module, only: nemsio_gfile,nemsio_open,nemsio_close,& - nemsio_getfilehead,nemsio_getheadvar,nemsio_realkind,nemsio_charkind,& - nemsio_readrecv,nemsio_init,nemsio_setheadvar,nemsio_writerecv - - use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& - quantize_data, read_attribute, close_dataset, get_dim, read_vardata - - implicit none - - character(len=500) :: filename - - real(r_kind), dimension(npts,nlevs+1) :: pressi -! real(r_kind), dimension(nlons*nlats) :: ug - real(r_single), dimension(npts) :: tmpgrd - ! type(sigio_data) sigdata - type(nemsio_gfile) :: gfile - real(nemsio_realkind), dimension(nlons*nlats) :: nems_wrk - real(r_kind), dimension(nlons*nlats) :: psfc - real(r_single),allocatable,dimension(:,:,:) :: nems_vcoord - real(r_kind), allocatable, dimension(:) :: ak,bk - real(r_kind) :: sumcoslat - - integer(i_kind) :: nlevsin,nlonsin,nlatsin,idvc - integer(i_kind) :: i,j,k,iret!iunitsig,iret - - ! newly added on 2020.07.09 - type(Dataset) :: dset - type(Dimension) :: londim,latdim,levdim - integer(i_kind) :: iunitsig - type(sigio_head) :: sighead - type(sigio_data) :: sigdata - real(r_kind), allocatable, dimension(:) :: psg - real(r_kind), dimension(ndimspec) :: vrtspec - real(r_single), allocatable, dimension(:,:) :: values_2d - integer(i_kind) :: ierr - - - - iunitsig = 77 - - ! ============================================================ - ! Read analysis data - ! ================================== - ! 2020.07.09 by LL - ! update ncio - !if (nproc .eq. 0) then - - filename = trim(adjustl(datapath))//trim(adjustl(andataname)) - if (nproc == 0) print *,'reading analysis file: ',filename - ! --- nemsio data ------------------------------------------------- - if (use_gfs_nemsio) then - call nemsio_init(iret=iret) - if(iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_init,iret=',iret - call stop2(23) - end if - - call nemsio_open(gfile,filename,'READ',iret=iret) - - if (iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_open,iret=',iret - call stop2(23) - endif - - call nemsio_getfilehead(gfile,iret=iret, dimx=nlonsin, dimy=nlatsin,& - dimz=nlevsin,idvc=idvc) - if (nlons /= nlonsin .or. nlats /= nlatsin .or. nlevs /= nlevsin) then - print *,'incorrect dims in nemsio file' - print *,'expected',nlons,nlats,nlevs - print *,'got',nlonsin,nlatsin,nlevsin - call stop2(23) - end if - - - ! --- NetCDF data ------------------------------------------------- - else if (use_gfs_ncio) then - - - dset = open_dataset(filename) - - londim = get_dim(dset,'grid_xt'); nlonsin = londim%len - latdim = get_dim(dset,'grid_yt'); nlatsin = latdim%len - levdim = get_dim(dset,'pfull'); nlevsin = levdim%len - idvc=2 - - - ! --- Other data ------------------------------------------------- - else - call sigio_srohdc(iunitsig,trim(filename), & - sighead,sigdata,iret) - if (iret /= 0) then - print *,'error reading file in gridio ',trim(filename) - call stop2(23) - end if - endif - - !endif - ! ============================================================ - - allocate(ak(nlevs+1)) - allocate(bk(nlevs+1)) - allocate(psg(nlons*nlats)) - allocate(weight(npts,nlevs)) - allocate(grweight(npts)) - - - if (.not. constants_initialized) then - print *,'constants not initialized (with init_constants, init_constants_derived)' - call stop2(23) - end if - - ! calculate weight on the grid point - sumcoslat = zero - - ! if(reducedgrid) then - k=0 - do i=1,nlats - do j=1,lonsperlat(i) - k=k+1 - grweight(k) = cos(gaulats(i)) * real(nlonsfull,r_kind) & - / real(lonsperlat(i),r_kind) - sumcoslat = sumcoslat + grweight(k) - end do - end do - - ! else - ! do i=1,nlons*nlats - ! grweight(i) = cos(latsgrd(i)) - ! sumcoslat = sumcoslat + grweight(i) - ! end do - ! end if - - sumcoslat = 1.0_r_kind / sumcoslat - grweight(:) = sqrt(grweight(:)*sumcoslat) - - ! ==================================================== - ! Extract surface pressure in pa - ! to aid in quantifying mass - ! ======================================== - ! === Option ONE, nemsio === - if (use_gfs_nemsio) then - - - call nemsio_readrecv(gfile,'pres','sfc',1,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_readrecv(ps), iret=',iret - call stop2(23) - endif - - ! Assign surface pressure in pa - psfc = nems_wrk - - ! Extract surface pressure - ! on reduced gaussian grid - call regtoreduced(psfc,tmpgrd) - - ! calculate half level pressures - if (allocated(nems_vcoord)) deallocate(nems_vcoord) - - allocate(nems_vcoord(nlevs+1,3,2)) - call nemsio_getfilehead(gfile,iret=iret,vcoord=nems_vcoord) - - if ( iret /= 0 ) then - write(6,*)' gridio: ***ERROR*** problem reading header ', & - 'vcoord, Status = ',iret - call stop2(99) - endif - - if ( idvc == 0 ) then ! sigma coordinate, old file format. - ak = zero - bk = nems_vcoord(1:nlevs+1,1,1) - else if ( idvc == 1 ) then ! sigma coordinate - ak = zero - bk = nems_vcoord(1:nlevs+1,2,1) - - else if ( idvc == 2 .or. idvc == 3 ) then ! hybrid coordinate - - ! AFE ak = nems_vcoord(1:nlevs+1,1,1) - ! AFE ak = nems_vcoord(1:nlevs+1,2,1) - ak = 0.01_r_kind*nems_vcoord(1:nlevs+1,1,1) ! convert to mb - bk = nems_vcoord(1:nlevs+1,2,1) - else - write(6,*)'gridio: ***ERROR*** INVALID value for idvc=',idvc - call stop2(85) - end if - - !==> pressure at interfaces. - do k=1,nlevs+1 - pressi(:,k)=ak(k)+bk(k)*tmpgrd - end do - deallocate(ak,bk) - - ! === Option TWO, NetCDF === - else if (use_gfs_ncio) then - - call mpi_barrier(mpi_comm_world,ierr) - - call read_vardata(dset, 'pressfc', values_2d,errcode=iret) - - if (iret /= 0) then - print *,'error reading ps' - call stop2(31) - endif - - psg = 0.01_r_kind*reshape(values_2d,(/nlons*nlats/)) - - ! Extract surface pressure - ! on reduced gaussian grid - call regtoreduced(psg,tmpgrd) - - call read_attribute(dset, 'ak', ak) - call read_attribute(dset, 'bk', bk) - - ! pressure at interfaces - do k=1,nlevs+1 - pressi(:,k) = 0.01_r_kind*ak(nlevs-k+2)+bk(nlevs-k+2)*tmpgrd - - if (nproc == 0) print *,'netcdf, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) - enddo - - deallocate(ak,bk,values_2d) - - ! === Option THREE, other === - else - - vrtspec = sigdata%ps - call sptez_s(vrtspec,psg,1) - !==> input psg is ln(ps) in centibars - convert to ps in millibars. - psg = 10._r_kind*exp(psg) - allocate(ak(nlevs+1),bk(nlevs+1)) - if (sighead%idvc .eq. 0) then ! sigma coordinate, old file format. - ak = zero - bk = sighead%si(1:nlevs+1) - else if (sighead%idvc == 1) then ! sigma coordinate - ak = zero - bk = sighead%vcoord(1:nlevs+1,2) - else if (sighead%idvc == 2 .or. sighead%idvc == 3) then ! hybrid coordinate - bk = sighead%vcoord(1:nlevs+1,2) - ak = 0.01_r_kind*sighead%vcoord(1:nlevs+1,1) ! convert to mb - else - print *,'unknown vertical coordinate type',sighead%idvc - call stop2(32) - end if - - do k=1,nlevs+1 - pressi(:,k)=ak(k)+bk(k)*psg - print *,'sigio, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) - enddo - - deallocate(ak,bk) - - endif - -!!! LL - if (use_gfs_nemsio) call nemsio_close(gfile,iret=iret) - if (use_gfs_ncio) call close_dataset(dset) - - - - - !$omp parallel do private(k) shared(weight,pressi,grweight,nlevs) - do k=1,nlevs - ! sqrt(dp)*sqrt(area) - weight(:,k)=sqrt( (pressi(:,k)-pressi(:,k+1))/tmpgrd(:))*grweight(:) - end do - - !$omp end parallel do - - - - - - - - return - - end subroutine get_weight - - - - - - - - - subroutine destroy_weight - implicit none - if(allocated(weight)) deallocate(weight) - if(allocated(grweight)) deallocate(grweight) - end subroutine destroy_weight - - - - - - subroutine divide_weight(grdin) - implicit none - real(r_single), dimension(npts_max,ncdim), intent(inout) :: grdin - real(r_single) cptr,qweight,rdtrpr - integer(i_kind) :: k,npt - cptr = real(sqrt(tref/cp),r_kind) - qweight = real(sqrt(cp*tref/wmoist)/hvap,r_kind) - rdtrpr = real(sqrt(pref/(rd*tref)),r_kind) - do npt=1,numptsperproc(nproc+1) - do k=1,nlevs - grdin(npt,k) = grdin(npt,k) / weight(indxproc(nproc+1,npt),k) - grdin(npt,nlevs+k) = grdin(npt,nlevs+k) / weight(indxproc(nproc+1,npt),k) - grdin(npt,2*nlevs+k) = grdin(npt,2*nlevs+k) * cptr / weight(indxproc(nproc+1,npt),k) - if (nvars .gt. 3) then - grdin(npt,3*nlevs+k) = grdin(npt,3*nlevs+k) * qweight / weight(indxproc(nproc+1,npt),k) - end if - end do -! AFE the indexing schema needs to be cleaned up -! grdin(npt,nvars*nlevs+1) = grdin(npt,nvars*nlevs+1) & - grdin(npt,ncdim) = grdin(npt,ncdim) & - & * rdtrpr / grweight(indxproc(nproc+1,npt)) - end do - return - end subroutine divide_weight - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine readgriddata_efsoi(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,mode,nanal,ft,hr,infilename) - - use sigio_module, only: sigio_head, sigio_data, sigio_sclose, sigio_sropen, & - sigio_srohdc, sigio_sclose, sigio_aldata, sigio_axdata - - use nemsio_module, only: nemsio_gfile,nemsio_open,nemsio_close,& - nemsio_getfilehead,nemsio_getheadvar,nemsio_realkind,nemsio_charkind,& - nemsio_readrecv,nemsio_init,nemsio_setheadvar,nemsio_writerecv - - use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& - read_attribute, close_dataset, get_dim, read_vardata - - implicit none - - integer, intent(in), optional :: nanal - integer, intent(in), optional :: ft - integer, intent(in), optional :: hr - -! character, intent(in), optional :: infilename AFE - character(len=100), intent(in), optional :: infilename ! AFE - integer, intent(in) :: mode - character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d - character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d - character(len=7) charnanal - integer, intent(in) :: n2d, n3d - integer, dimension(0:n3d), intent(in) :: levels - integer, intent(in) :: ndim - real(r_single), dimension(npts,ndim), intent(out) :: grdin - real(r_kind) cptr,qweight,rdtrpr - character(len=500) :: filename - character(len=3) charft - character(len=2) charhr - - real(r_kind), dimension(nlons*nlats) :: ug,vg - real(r_single), dimension(npts,nlevs) :: tv, q, tv_to_t - real(r_kind), allocatable, dimension(:) :: psg - - real(r_single),allocatable,dimension(:,:,:) :: nems_vcoord - real(nemsio_realkind), dimension(nlons*nlats) :: nems_wrk,nems_wrk2 - type(nemsio_gfile) :: gfile - - integer(i_kind) :: u_ind, v_ind, tv_ind, q_ind - integer(i_kind) :: ps_ind - - integer(i_kind) :: k,iret,idvc,nlonsin,nlatsin,nlevsin - - ! Added by LL - character(len=10) :: fileformat - integer(i_kind) :: iunitsig - type(sigio_head) :: sighead - type(sigio_data) :: sigdata - type(Dataset) :: dset - type(Dimension) :: londim,latdim,levdim - real(r_kind), dimension(ndimspec) :: vrtspec,divspec - real(r_single), allocatable, dimension(:,:,:) :: ug3d,vg3d - real(r_single), allocatable, dimension(:,:) :: values_2d - integer(i_kind) :: nb,ne - - - - ! Added by LL - iunitsig = 77 - nb = 1 - ne = 1 - - ! ----------------------------! - ! EFSOI Filename constructions! - ! ----------------------------! - - ! Construct the Filename based on - ! input arguments - if(.not. present(infilename)) then - write(charft, '(i3.3)') ft - write(charhr, '(i2.2)') hr - if (nanal > 0) then - write(charnanal,'(a3, i3.3)') 'mem', nanal - else - charnanal = 'ensmean' - endif - endif - - ! ====================================== - if (use_gfs_nemsio) then - fileformat = '.nemsio' - else if (use_gfs_ncio) then - fileformat = '.nc' - else - print *,'Warning in gridio_efsoi.f90 by LL (2020.07.09)' - end if - ! ====================================== - - - ! === EFSOI filename construction ============= - if(forecast_impact) then - - if(mode == read_ensmean_forecast) then - filename = trim(adjustl(datapath))//"gdas.t"//charhr//"z.atmf"//charft// & - "."//trim(adjustl(charnanal))//trim(fileformat) - else if(mode == read_analysis_mean) then - filename = trim(adjustl(datapath))//"gdas.t"//charhr//"z.atmanl."// & - trim(adjustl(charnanal))//trim(fileformat) - else if(mode == read_member_forecasts) then - filename = trim(adjustl(datapath))//trim(adjustl(charnanal))//"/"// & - "gdas.t"//charhr//"z.atmf"//charft//trim(fileformat) - else if(mode == read_verification) then - filename = trim(adjustl(datapath))//infilename - else - print *,'This mode is not supported: mode=',mode - call stop2(23) - end if - else - ! Analysis Impact - if(mode == read_ensmean_forecast) then - filename = trim(adjustl(datapath))//"gdas.t"//charhr//"z.atmf"//charft// & - "."//trim(adjustl(charnanal))//trim(fileformat) - else if(mode == read_member_analyses) then - filename = trim(adjustl(datapath))//trim(adjustl(charnanal))//"/"//"gdas.t"// & - charhr//"z.atmanl"//trim(fileformat) - else if(mode == read_verification) then - filename = trim(adjustl(datapath))//infilename - else - print *,'This mode is not supported: mode=',mode - call stop2(23) - end if - endif - - ! --------------------------------------------! - ! Read in state vector from file and transform! - ! to EFSOI relevant quantities ! - ! --------------------------------------------! - ! Updated by LL on 2020.07.09 - - - if (nproc == 0) print *,'reading state vector file: ',filename - - if (use_gfs_nemsio) then - - call nemsio_init(iret=iret) - if(iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_init, iret=',iret - call stop2(23) - end if - - call nemsio_open(gfile,filename,'READ',iret=iret) - - if (iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_open, iret=',iret - call stop2(23) - endif - call nemsio_getfilehead(gfile,iret=iret, dimx=nlonsin, dimy=nlatsin,& - dimz=nlevsin,idvc=idvc) - if (nlons /= nlonsin .or. nlats /= nlatsin .or. nlevs /= nlevsin) then - print *,'incorrect dims in nemsio file' - print *,'expected',nlons,nlats,nlevs - print *,'got',nlonsin,nlatsin,nlevsin - call stop2(23) - end if - - else if (use_gfs_ncio) then - - ! Edited by LL on 2020.07.16 - dset = open_dataset(filename) - - londim = get_dim(dset,'grid_xt'); nlonsin = londim%len - latdim = get_dim(dset,'grid_yt'); nlatsin = latdim%len - levdim = get_dim(dset,'pfull'); nlevsin = levdim%len - idvc=2 - - else - - call sigio_srohdc(iunitsig,trim(filename), & - sighead,sigdata,iret) - if (iret /= 0) then - print *,'error reading file in gridio ',trim(filename) - call stop2(23) - end if - endif - - - cptr = sqrt(cp/tref) - qweight = sqrt(wmoist/(cp*tref))*hvap - rdtrpr = sqrt(rd*tref)/pref - - - u_ind = getindex(vars3d, 'u') !< indices in the state var arrays - v_ind = getindex(vars3d, 'v') ! U and V (3D) - tv_ind = getindex(vars3d, 'tv') ! Tv (3D) - q_ind = getindex(vars3d, 'q') ! Q (3D) - ps_ind = getindex(vars2d, 'ps') ! Ps (2D) - - - - if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4) - - allocate(psg(nlons*nlats)) - - - ! ====================================================== - ! Get surface pressure - ! ==================== - if (use_gfs_nemsio) then - - call nemsio_readrecv(gfile,'pres','sfc',1,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_readrecv(ps), iret=',iret - call stop2(23) - endif - psg = 0.01_r_kind*nems_wrk ! convert ps to millibars. - - ! - if (allocated(nems_vcoord)) deallocate(nems_vcoord) - allocate(nems_vcoord(nlevs+1,3,2)) - call nemsio_getfilehead(gfile,iret=iret,vcoord=nems_vcoord) - if ( iret /= 0 ) then - write(6,*)' gridio: ***ERROR*** problem reading header ', & - 'vcoord, Status = ',iret - call stop2(99) - endif - - else if (use_gfs_ncio) then - - call read_vardata(dset, 'pressfc', values_2d,errcode=iret) - if (iret /= 0) then - print *,'error reading ps' - call stop2(31) - endif - psg = 0.01_r_kind*reshape(values_2d,(/nlons*nlats/)) - - else - vrtspec = sigdata%ps - call sptez_s(vrtspec,psg,1) - !==> input psg is ln(ps) in centibars - convert to ps in millibars. - psg = 10._r_kind*exp(psg) - endif - - - - ! ============================================================================= - ! get U,V,temp,q,ps on gaussian grid. - ! ===================================== - ! u is first nlevs, v is second, t is third, then tracers. - if (use_gfs_nemsio) then - - do k=1,nlevs - - ! Get u-wind - call nemsio_readrecv(gfile,'ugrd','mid layer',k,nems_wrk,iret=iret) - - if (iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_readrecv(ugrd), iret=',iret - call stop2(23) - endif - - ug = nems_wrk - - ! Get v-wind - call nemsio_readrecv(gfile,'vgrd','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_readrecv(vgrd), iret=',iret - call stop2(23) - endif - - vg = nems_wrk - - if (u_ind > 0) call copytogrdin(ug,grdin(:,levels(u_ind-1) + k)) - if (v_ind > 0) call copytogrdin(vg,grdin(:,levels(v_ind-1) + k)) - - - ! Transformation to EFSOI relevant quantities - ! Assign weighted kinetic energy components. There - ! are no unit/metric differences for the kinetic component - grdin(:,levels(u_ind-1) + k) = weight(:,k) * grdin(:,levels(u_ind-1) + k) - grdin(:,levels(v_ind-1) + k) = weight(:,k) * grdin(:,levels(v_ind-1) + k) - - call nemsio_readrecv(gfile,'tmp','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_readrecv(tmp), iret=',iret - call stop2(23) - endif - - call nemsio_readrecv(gfile,'spfh','mid layer',k,nems_wrk2,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata_efsoi: GFS: problem with nemsio_readrecv(spfh), iret=',iret - call stop2(23) - endif - - ug = nems_wrk - vg = nems_wrk2 - - call copytogrdin(ug,tv(:,k)) - call copytogrdin(vg, q(:,k)) - - ! Transformation to EFSOI relevant quantities - ! Mass component quantities - tv(:,k) = cptr * weight(:,k) * tv(:,k) - - ! Moisture component transormation to EFSOI - ! relevant quantities - q(:,k) = qweight * weight(:,k) * q(:,k) - tv_to_t(:,k) = ( one / (one + q(:,k)*0.61_r_kind) ) - - ! Approximate the necessary transformation - ! of virtual temperature to temperature - tv(:,k) = tv(:,k) * tv_to_t(:,k) - - if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k) = tv(:,k) - if (q_ind > 0) grdin(:,levels( q_ind-1)+k) = q(:,k) - - enddo - - else if (use_gfs_ncio) then - - ! ==== Get U and V ======================================== - call read_vardata(dset, 'ugrd', ug3d,errcode=iret) - - if (iret /= 0) then - print *,'error reading ugrd' - call stop2(22) - endif - - call read_vardata(dset, 'vgrd', vg3d,errcode=iret) - if (iret /= 0) then - print *,'error reading vgrd' - call stop2(23) - endif - - - do k=1,nlevs - ug = reshape(ug3d(:,:,nlevs-k+1),(/nlons*nlats/)) - vg = reshape(vg3d(:,:,nlevs-k+1),(/nlons*nlats/)) - - if (u_ind > 0) call copytogrdin(ug,grdin(:,levels(u_ind-1) + k)) - if (v_ind > 0) call copytogrdin(vg,grdin(:,levels(v_ind-1) + k)) - - ! Added by LL - ! Transformation to EFSOI relevant quantities - ! Assign weighted kinetic energy components. There - ! are no unit/metric differences for the kinetic component - grdin(:,levels(u_ind-1) + k) = weight(:,k) * grdin(:,levels(u_ind-1) + k) - grdin(:,levels(v_ind-1) + k) = weight(:,k) * grdin(:,levels(v_ind-1) + k) - - ! calculate vertical integral of mass flux div (ps tendency) - ! this variable is analyzed in order to enforce mass balance in the analysis - !if (pst_ind > 0) then - ! ug = ug*(pressi(:,k)-pressi(:,k+1)) - ! vg = vg*(pressi(:,k)-pressi(:,k+1)) - ! call sptezv_s(divspec,vrtspec,ug,vg,-1) ! u,v to div,vrt - ! call sptez_s(divspec,vmassdiv(:,k),1) ! divspec to divgrd - !endif - enddo - - - - ! ==== Get T and Q =========================================== - call read_vardata(dset,'tmp', ug3d,errcode=iret) - - if (iret /= 0) then - print *,'error reading tmp' - call stop2(24) - endif - - call read_vardata(dset,'spfh', vg3d,errcode=iret) - - if (iret /= 0) then - print *,'error reading spfh' - call stop2(25) - endif - - do k=1,nlevs - ug = reshape(ug3d(:,:,nlevs-k+1),(/nlons*nlats/)) - vg = reshape(vg3d(:,:,nlevs-k+1),(/nlons*nlats/)) - - !if (tsen_ind > 0) call copytogrdin(ug,grdin(:,levels(tsen_ind-1)+k,nb,ne)) - - call copytogrdin(vg, q(:,k)) - - ug = ug * ( 1.0 + fv*vg ) ! convert T to Tv - ! - call copytogrdin(ug,tv(:,k)) - - ! Transformation to EFSOI relevant quantities - ! Mass component quantities - tv(:,k) = cptr * weight(:,k) * tv(:,k) - ! Moisture component transormation to EFSOI - ! relevant quantities - q(:,k) = qweight * weight(:,k) * q(:,k) - tv_to_t(:,k) = ( one / (one + q(:,k)*0.61_r_kind) ) - - ! Approximate the necessary transformation - ! of virtual temperature to temperature - tv(:,k) = tv(:,k) * tv_to_t(:,k) - - if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k) = tv(:,k) - if (q_ind > 0) grdin(:,levels( q_ind-1)+k) = q(:,k) - enddo - - deallocate(ug3d,vg3d) - - else - - ! $omp parallel do private(k,ug,vg,divspec,vrtspec) - ! shared(sigdata,pressi,vmassdiv,grdin,tv,q,cw,u_ind,v_ind,pst_ind,q_ind,tsen_ind,cw_ind,qi_ind,ql_ind) - do k=1,nlevs - - vrtspec = sigdata%z(:,k); divspec = sigdata%d(:,k) - call sptezv_s(divspec,vrtspec,ug,vg,1) - if (u_ind > 0) then - call copytogrdin(ug,grdin(:,levels(u_ind-1)+k)) - endif - if (v_ind > 0) then - call copytogrdin(vg,grdin(:,levels(v_ind-1)+k)) - endif - - ! calculate vertical integral of mass flux div (ps tendency) - ! this variable is analyzed in order to enforce mass balance in the analysis - !if (pst_ind > 0) then - ! ug = ug*(pressi(:,k)-pressi(:,k+1)) - ! vg = vg*(pressi(:,k)-pressi(:,k+1)) - ! call sptezv_s(divspec,vrtspec,ug,vg,-1) ! u,v to div,vrt - ! call sptez_s(divspec,vmassdiv(:,k),1) ! divspec to divgrd - !endif - - divspec = sigdata%t(:,k) - call sptez_s(divspec,ug,1) - call copytogrdin(ug,tv(:,k)) - - divspec = sigdata%q(:,k,1) - call sptez_s(divspec,vg,1) - call copytogrdin(vg,q(:,k)) - - ! Transformation to EFSOI relevant quantities - ! Mass component quantities - tv(:,k) = cptr * weight(:,k) * tv(:,k) - ! Moisture component transormation to EFSOI - ! relevant quantities - q(:,k) = qweight * weight(:,k) * q(:,k) - tv_to_t(:,k) = ( one / (one + q(:,k)*0.61_r_kind) ) - - ! Approximate the necessary transformation - ! of virtual temperature to temperature - tv(:,k) = tv(:,k) * tv_to_t(:,k) - - if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k) = tv(:,k) - if (q_ind > 0) grdin(:,levels( q_ind-1)+k) = q(:,k) - - enddo - endif - - - ! ============================================================= - ! surface pressure - ! ================ - if (ps_ind > 0) then - call copytogrdin(psg,grdin(:,levels(n3d) + ps_ind)) - endif - - ! Transformation to EFSOI relevant quantities - ! Surface pressure contribution - grdin(:,levels(n3d) + ps_ind) = rdtrpr * grweight(:) * 100._r_kind * grdin(:,levels(n3d) + ps_ind) - - - - - deallocate(psg) - - ! updated by LL - if (use_gfs_nemsio) call nemsio_close(gfile,iret=iret) - if (use_gfs_ncio) call close_dataset(dset) - - return - - contains - - ! copying to grdin (calling regtoreduced if reduced grid) - subroutine copytogrdin(field, grdin) - implicit none - - real(r_kind), dimension(:), intent(in) :: field - real(r_single), dimension(:), intent(inout) :: grdin - - call regtoreduced(field, grdin) - - end subroutine copytogrdin - - end subroutine readgriddata_efsoi - -end module gridio_efsoi diff --git a/util/EFSOI_Utilities/src/loadbal_efsoi.f90 b/util/EFSOI_Utilities/src/loadbal_efsoi.f90 deleted file mode 100644 index 1a9f046a64..0000000000 --- a/util/EFSOI_Utilities/src/loadbal_efsoi.f90 +++ /dev/null @@ -1,381 +0,0 @@ -module loadbal_efsoi -!$$$ module documentation block -! -! module: loadbal decompose ob priors and horizontal grid points -! to minimize load imbalance. Creates various -! arrays that define the decomposition. -! -! prgmmr: whitaker org: esrl/psd date: 2009-02-23 -! groff org: emc date: 2018-09-03 -! abstract: -! -! Public Subroutines: -! load_balance_efsoi: set up decomposition (for ob priors and analysis grid points) -! that minimizes load imbalance. -! The decomposition uses "Graham's rule", which simply -! stated, assigns each new work item to the task that currently has the -! smallest load. -! loadbal_cleanup_efsoi: deallocate allocated arrays. -! -! Private Subroutines: -! estimate_work_efsoi: estimate work needed to update each analysis grid -! point (considering all the observations within the localization radius). -! Public Variables (all defined by subroutine load_balance_efsoi): -! npts_min: (integer scalar) smallest number of grid points assigned to a task. -! npts_max: (integer scalar) maximum number of grid points assigned to a task. -! nobs_min: (integer scalar, serial enkf only) smallest number of observation priors assigned to a task. -! nobs_max: (integer scalar, serial enkf only) maximum number of observation priors assigned to a task. -! numproc: (integer scalar) total number of MPI tasks (from module mpisetup) -! nobstot: (integer scalar) total number of obs to be assimilated (from module -! enkf_obsmod). -! numobsperproc(numproc): (serial enkf only) integer array containing # of ob priors -! assigned to each task. -! numptsperproc(numproc): integer array containing # of grid points assigned to -! each task. -! indxproc(numproc,npts_max): integer array with the indices (1,2,...npts) of -! analysis grid points assigned to each task. -! indxproc_obs(numproc,nobs_max): (serial enkf only) integer array with the indices -! (1,2,...nobstot) of observation priors assigned to that task. -! iprocob(nobstot): (serial enkf only) integer array containing the task number that has been -! assigned to update each observation prior. -! indxob_chunk(nobstot): (serial enkf only) integer array that maps the index of the ob priors -! being assimilated (1,2,3...nobstot) to the index of the obs being -! updated on that task (1,2,3,...numobsperproc(nproc)) - inverse of -! indxproc_obs. -! ensmean_obchunk(nobs_max): (serial enkf only) real array of ensemble mean observation priors -! being updated on that task (use indxproc_obs to find -! corresponding index in ensemble_ob). -! obloc_chunk(3,nobs_max): (serial enkf only) real array of spherical cartesian coordinates -! of ob priors being updated on this task. -! grdloc_chunk(3,npts_max): real array of spherical cartesian coordinates -! of analysis grid points being updated on this task. -! lnp_chunk(npts_max,ncdim): real array of log(pressures) of control variables -! being updated on this task. -! oblnp_chunk(nobs_max,ncdim): (serial enkf only) real array of log(pressures) of ob priors -! being updated on this task. -! obtime_chunk(nobs_max): (serial enkf only) real array of ob times of ob priors -! being updated on this task (expressed as an offset from the analysis time in -! hours). -! anal_obchunk_prior(nanals,nobs_max): (serial enkf only) real array of observation prior -! ensemble perturbations to be updated on this task (not used in LETKF). -! kdtree_grid: pointer to kd-tree structure used for nearest neighbor searches -! for model grid points (only searches grid points assigned to this task). -! kdtree_obs: pointer to kd-tree structure used for nearest neighbor searches -! for observations (only searches ob locations assigned to this task). -! kdtree_obs2: (LETKF only) pointer to kd-tree structure used for nearest neighbor searches -! for observations (searches all observations) -! anal_chunk(nanals,npts_max,ncdim,nbackgrounds): real array of ensemble perturbations -! updated on each task. -! anal_chunk_prior(nanals,npts_max,ncdim,nbackgrounds): real array of prior ensemble -! perturbations. Before analysis anal_chunk=anal_chunk_prior, after -! analysis anal_chunk contains posterior perturbations. -! ensmean_chunk(npts_max,ncdim,nbackgrounds): real array containing pieces of ensemble -! mean to be updated on each task. -! ensmean_chunk_prior(npts_max,ncdim,nbackgrounds): as above, for ensemble mean prior. -! Before analysis ensmean_chunk=ensmean_chunk_prior, after analysis -! ensmean_chunk contains posterior ensemble mean. -! -! -! Modules Used: mpisetup, params, kinds, constants, enkf_obsmod, gridinfo, -! kdtree_module, covlocal, controlvec -! -! program history log: -! 2009-02-23 Initial version. -! 2011-06-21 Added the option of observation box selection for LETKF. -! 2015-07-25 Remove observation box selection (use kdtree instead). -! 2016-05-02 shlyaeva: modification for reading state vector from table -! 2016-09-07 shlyaeva: moved distribution of chunks here from controlvec -! -! attributes: -! language: f95 -! -!$$$ - -use mpisetup -use params, only: datapath, nanals, simple_partition, & - corrlengthnh, corrlengthsh, corrlengthtr, lupd_obspace_serial,& - efsoi_flag -use enkf_obsmod, only: nobstot, obloc, oblnp, ensmean_ob, obtime, anal_ob, corrlengthsq -use kinds, only: r_kind, i_kind, r_double, r_single -use kdtree2_module, only: kdtree2, kdtree2_create, kdtree2_destroy, & - kdtree2_result, kdtree2_r_nearest -use gridinfo_efsoi, only: gridloc, logp, latsgrd, nlevs_pres, npts -use constants, only: zero, rad2deg, deg2rad - -implicit none -private -public :: load_balance_efsoi, loadbal_cleanup_efsoi - -real(r_single),public, allocatable, dimension(:,:) :: lnp_chunk, & - anal_obchunk_prior -real(r_single),public, allocatable, dimension(:,:,:) :: anal_chunk, anal_chunk_prior -real(r_single),public, allocatable, dimension(:,:) :: ensmean_chunk, ensmean_chunk_prior - -! arrays passed to kdtree2 routines need to be single -real(r_single),public, allocatable, dimension(:,:) :: obloc_chunk, grdloc_chunk -real(r_single),public, allocatable, dimension(:) :: oblnp_chunk, & - obtime_chunk, ensmean_obchunk -integer(i_kind),public, allocatable, dimension(:) :: iprocob, indxob_chunk,& - numptsperproc, numobsperproc -integer(i_kind),public, allocatable, dimension(:,:) :: indxproc, indxproc_obs -integer(i_kind),public :: npts_min, npts_max, nobs_min, nobs_max -integer(8) totsize -! kd-tree structures. -type(kdtree2),public,pointer :: kdtree_obs, kdtree_grid, kdtree_obs2 - -contains - -subroutine load_balance_efsoi() -! set up decomposition (for analysis grid points, and ob priors in serial EnKF) -! that minimizes load imbalance. -! Uses "Graham's rule", which simply -! stated, assigns each new work item to the task that currently has the -! smallest load. -implicit none -integer(i_kind), allocatable, dimension(:) :: rtmp,numobs -!real(r_single), allocatable, dimension(:) :: buffer -integer(i_kind) np,i,n,nn,nob1,nob2,ierr -real(r_double) t1 -logical test_loadbal - -! partition state vector for using Grahams rule.. -! ("When a new job arrives, allocate it to the server -! that currently has the smallest load") -allocate(numobs(npts)) -allocate(numptsperproc(numproc)) -allocate(rtmp(numproc)) -t1 = mpi_wtime() -! assume work load proportional to number of 'nearby' obs -call estimate_work_efsoi(numobs) ! fill numobs array with number of obs per horiz point -! distribute the results of estimate_work to all processors. -call mpi_allreduce(mpi_in_place,numobs,npts,mpi_integer,mpi_sum,mpi_comm_world,ierr) -if (nproc == 0) print *,'time in estimate_work_efsoi = ',mpi_wtime()-t1,' secs' -if (nproc == 0) print *,'min/max numobs',minval(numobs),maxval(numobs) -! loop over horizontal grid points on analysis grid. -t1 = mpi_wtime() -rtmp = 0 -numptsperproc = 0 -np = 0 -test_loadbal = .false. ! simple partition for testing -! -!PRINT *, npts, 'npts load bal' -do n=1,npts - if (test_loadbal) then - ! use simple partition (does not use estimated workload) for testing - np = np + 1 - if (np > numproc) np = 1 - else - np = minloc(rtmp,dim=1) - ! np is processor with the fewest number of obs to process - ! add this grid point to list for nmin - rtmp(np) = rtmp(np)+numobs(n) - endif - numptsperproc(np) = numptsperproc(np)+1 - !PRINT *, numptsperproc, 'loadbal nptsperpr' - - - -end do - - -npts_max = maxval(numptsperproc) -npts_min = minval(numptsperproc) - -allocate(indxproc(numproc,npts_max)) -! indxproc(np,i) is i'th horiz grid index for processor np. -! there are numptsperpoc(np) i values for processor np -rtmp = 0 -numptsperproc = 0 -np = 0 -do n=1,npts - if (test_loadbal) then - ! use simple partition (does not use estimated workload) for testing - np = np + 1 - if (np > numproc) np = 1 - else - np = minloc(rtmp,dim=1) - rtmp(np) = rtmp(np)+numobs(n) - endif - numptsperproc(np) = numptsperproc(np)+1 ! recalculate - indxproc(np,numptsperproc(np)) = n -end do - -!PRINT *, numptsperproc, 'second time' -! print estimated workload for each task -if (nproc == 0) then - do np=1,numproc - rtmp(np) = 0 - do n=1,numptsperproc(np) - rtmp(np) = rtmp(np) + numobs(indxproc(np,n)) - enddo - - - enddo - - print *,'min/max estimated work ',& - minval(rtmp),maxval(rtmp) - -endif -deallocate(rtmp,numobs) - -if (nproc == 0) then - print *,'npts = ',npts - print *,'min/max number of points per proc = ',npts_min,npts_max - print *,'time to do model space decomp = ',mpi_wtime()-t1 -end if - -! setup arrays to hold subsets of grid information for each task. -allocate(grdloc_chunk(3,numptsperproc(nproc+1))) -allocate(lnp_chunk(numptsperproc(nproc+1),nlevs_pres)) -do i=1,numptsperproc(nproc+1) - grdloc_chunk(:,i) = gridloc(:,indxproc(nproc+1,i)) - do nn=1,nlevs_pres - lnp_chunk(i,nn) = logp(indxproc(nproc+1,i),nn) - end do -end do - -allocate(numobsperproc(numproc)) -allocate(iprocob(nobstot)) -! default is to partition obs simply, since -! speed up from using Graham's rule for observation process -! often does not justify cost of estimating workload in ob space. - -! distribute obs without trying to estimate workload -t1 = mpi_wtime() -numobsperproc = 0 -np=0 -do n=1,nobstot - np=np+1 - if(np > numproc)np = 1 - numobsperproc(np) = numobsperproc(np)+1 - iprocob(n) = np-1 -enddo -nobs_min = minval(numobsperproc) -nobs_max = maxval(numobsperproc) -allocate(indxproc_obs(numproc,nobs_max)) -numobsperproc = 0 -do n=1,nobstot - np=iprocob(n)+1 - numobsperproc(np) = numobsperproc(np)+1 ! recalculate - ! indxproc_obs(np,i) is i'th ob index for processor np. - ! there are numobsperpoc(np) i values for processor np - indxproc_obs(np,numobsperproc(np)) = n -end do -if (nproc == 0) then - print *,'nobstot = ',nobstot - print *,'min/max number of obs per proc = ',nobs_min,nobs_max - print *,'time to do ob space decomp = ',mpi_wtime()-t1 -end if -! for serial enkf, send out observation priors to be updated on each processor. -allocate(anal_obchunk_prior(nanals,nobs_max)) -if(nproc == 0) then - print *,'sending out observation prior ensemble perts from root ...' - totsize = nobstot - totsize = totsize*nanals - print *,'nobstot*nanals',totsize - t1 = mpi_wtime() - ! send one big message to each task. - do np=1,numproc-1 - do nob1=1,numobsperproc(np+1) - nob2 = indxproc_obs(np+1,nob1) - anal_obchunk_prior(1:nanals,nob1) = anal_ob(1:nanals,nob2) - end do - call mpi_send(anal_obchunk_prior,nobs_max*nanals,mpi_real4,np, & - 1,mpi_comm_world,ierr) - end do - ! anal_obchunk_prior on root (no send necessary) - do nob1=1,numobsperproc(1) - nob2 = indxproc_obs(1,nob1) - anal_obchunk_prior(1:nanals,nob1) = anal_ob(1:nanals,nob2) - end do - ! now we don't need anal_ob anymore for serial EnKF. - if (.not. lupd_obspace_serial) deallocate(anal_ob) -else - ! recv one large message on each task. - call mpi_recv(anal_obchunk_prior,nobs_max*nanals,mpi_real4,0, & - 1,mpi_comm_world,mpi_status,ierr) -end if -call mpi_barrier(mpi_comm_world, ierr) -if(nproc == 0) print *,'... took ',mpi_wtime()-t1,' secs' -! these arrays only needed for serial filter -! nob1 is the index of the obs to be processed on this rank -! nob2 maps nob1 to 1:nobstot array (nobx) -allocate(obloc_chunk(3,numobsperproc(nproc+1))) -allocate(oblnp_chunk(numobsperproc(nproc+1))) -allocate(obtime_chunk(numobsperproc(nproc+1))) -allocate(ensmean_obchunk(numobsperproc(nproc+1))) -allocate(indxob_chunk(nobstot)) -indxob_chunk = -1 -do nob1=1,numobsperproc(nproc+1) - nob2 = indxproc_obs(nproc+1,nob1) - oblnp_chunk(nob1) = oblnp(nob2) - obtime_chunk(nob1) = obtime(nob2) - indxob_chunk(nob2) = nob1 - ensmean_obchunk(nob1) = ensmean_ob(nob2) - obloc_chunk(:,nob1) = obloc(:,nob2) -enddo -! set up kd-tree for efsoi to search only the subset -! of gridpoints, obs to be updated on this processor.. -if (numptsperproc(nproc+1) >= 3 .and. .not. lupd_obspace_serial .and. .not. efsoi_flag) then - kdtree_grid => kdtree2_create(grdloc_chunk,sort=.false.,rearrange=.true.) -endif -if (numobsperproc(nproc+1) >= 3) then - kdtree_obs => kdtree2_create(obloc_chunk,sort=.false.,rearrange=.true.) -end if - -end subroutine load_balance_efsoi - -subroutine estimate_work_efsoi(numobs) -! estimate work needed to update each analysis grid -! point (considering all the observations within the localization radius). -use covlocal, only: latval - -implicit none -integer(i_kind), dimension(:), intent(inout) :: numobs -real(r_single) :: deglat,corrlength,corrsq -type(kdtree2_result),dimension(:),allocatable :: sresults - -integer nob,n1,n2,i,ideln - -ideln = int(real(npts)/real(numproc)) -n1 = 1 + nproc*ideln -n2 = (nproc+1)*ideln -if (nproc == numproc-1) n2 = npts - -! loop over 'good' obs. -numobs = 1 ! set min # of obs to 1, not 0 (so single ob test behaves) -!$omp parallel do schedule(dynamic,1) private(nob,i,deglat,corrlength,sresults,corrsq) -obsloop: do i=n1,n2 - do nob=1,nobstot - if (sum((obloc(1:3,nob)-gridloc(1:3,i))**2,1) < corrlengthsq(nob)) & - numobs(i) = numobs(i) + 1 - end do -end do obsloop -!$omp end parallel do - -end subroutine estimate_work_efsoi - -subroutine loadbal_cleanup_efsoi() -! deallocate module-level allocatable arrays -if (allocated(anal_chunk)) deallocate(anal_chunk) -if (allocated(anal_chunk_prior)) deallocate(anal_chunk_prior) -if (allocated(ensmean_chunk)) deallocate(ensmean_chunk) -if (allocated(ensmean_chunk_prior)) deallocate(ensmean_chunk_prior) -if (allocated(obloc_chunk)) deallocate(obloc_chunk) -if (allocated(grdloc_chunk)) deallocate(grdloc_chunk) -if (allocated(lnp_chunk)) deallocate(lnp_chunk) -if (allocated(oblnp_chunk)) deallocate(oblnp_chunk) -if (allocated(obtime_chunk)) deallocate(obtime_chunk) -if (allocated(ensmean_obchunk)) deallocate(ensmean_obchunk) -if (allocated(iprocob)) deallocate(iprocob) -if (allocated(indxob_chunk)) deallocate(indxob_chunk) -if (allocated(indxproc_obs)) deallocate(indxproc_obs) -if (allocated(numptsperproc)) deallocate(numptsperproc) -if (allocated(numobsperproc)) deallocate(numobsperproc) -if (allocated(indxproc)) deallocate(indxproc) -if (associated(kdtree_obs)) call kdtree2_destroy(kdtree_obs) -if (associated(kdtree_obs2)) call kdtree2_destroy(kdtree_obs2) -if (associated(kdtree_grid)) call kdtree2_destroy(kdtree_grid) -end subroutine loadbal_cleanup_efsoi - -end module loadbal_efsoi diff --git a/util/EFSOI_Utilities/src/loc_advection.f90 b/util/EFSOI_Utilities/src/loc_advection.f90 deleted file mode 100644 index b6a544666d..0000000000 --- a/util/EFSOI_Utilities/src/loc_advection.f90 +++ /dev/null @@ -1,99 +0,0 @@ -module loc_advection - -use mpisetup -use gridinfo_efsoi, only: latsgrd,lonsgrd,nlevs_pres,npts,id_u,id_v -use loadbal_efsoi, only: numptsperproc, indxproc, kdtree_grid -use scatter_chunks_efsoi -use kinds -use params -use constants -use kdtree2_module -implicit none - -private -public loc_advection_efsoi, adloc_chunk -real(r_single),allocatable,dimension(:,:) :: adloc_chunk - -contains - -subroutine loc_advection_efsoi -!$$$ subprogram documentation block -! . . . . -! subprogram: loc_advection -! prgmmr: ota -! -! abstract: computes analysis grid point location with simple horizontal -! advection. -! -! program history log: -! 2011-09-16 ota - created -! -! input argument list: -! -! output argument list: -! -! attributes: -! language: f90 -! machine: -! -!$$$ end documentation block - implicit none - real(r_kind),allocatable,dimension(:) :: coslat - real(r_kind) :: adtime, halfpi, pi2, adlon, adlat - integer(i_kind) :: npt, nn, nnu, nnv, nnpt - allocate(adloc_chunk(3,numptsperproc(nproc+1)*nlevs_pres)) - allocate(coslat(numptsperproc(nproc+1))) - adtime = adrate*real(eft,r_kind)*3600._r_kind/rearth - halfpi = half * pi - pi2 = 2._r_kind * pi -!$omp parallel private(npt) -!$omp do - do npt=1,numptsperproc(nproc+1) - coslat(npt) = 1._r_kind/cos(latsgrd(indxproc(nproc+1,npt))) - end do -!$omp end do -!$omp end parallel -!$omp parallel private(nn,nnu,nnv,npt,nnpt,adlon,adlat) -!$omp do - do nn=1,nlevs_pres - if(nn == nlevs_pres) then - nnu = id_u(1) - nnv = id_v(1) - else - nnu = id_u(nn) - nnv = id_v(nn) - end if - do npt=1,numptsperproc(nproc+1) - nnpt = nlevs_pres * (npt-1) + nn - adlon = lonsgrd(indxproc(nproc+1,npt)) & - & - half * (ensmean_chunk(npt,nnu) + analmean_chunk(npt,nnu)) & - & * coslat(npt) * adtime - adlat = latsgrd(indxproc(nproc+1,npt)) & - & - half * (ensmean_chunk(npt,nnv) + analmean_chunk(npt,nnv)) & - & * adtime - if(adlat > halfpi) then - adlat = pi - adlat - adlon = adlon + pi - else if(adlat < -halfpi) then - adlat = -pi - adlat - adlon = adlon + pi - end if - if(adlon > pi2) then - adlon = mod(adlon,pi2) - else if(adlon < zero) then - adlon = mod(adlon,pi2) + pi2 - end if - adloc_chunk(1,nnpt) = cos(adlat)*cos(adlon) - adloc_chunk(2,nnpt) = cos(adlat)*sin(adlon) - adloc_chunk(3,nnpt) = sin(adlat) - end do - end do -!$omp end do -!$omp end parallel - deallocate(coslat) - ! kd-tree - kdtree_grid => kdtree2_create(adloc_chunk,sort=.false.,rearrange=.true.) - return -end subroutine loc_advection_efsoi - -end module loc_advection diff --git a/util/EFSOI_Utilities/src/scatter_chunks_efsoi.f90 b/util/EFSOI_Utilities/src/scatter_chunks_efsoi.f90 deleted file mode 100644 index 86e32ac41a..0000000000 --- a/util/EFSOI_Utilities/src/scatter_chunks_efsoi.f90 +++ /dev/null @@ -1,159 +0,0 @@ -module scatter_chunks_efsoi - -use mpisetup !, only: numproc, nproc, mpi_real4 -use kinds, only: i_kind, r_single -use loadbal_efsoi, only: npts_max, numptsperproc, indxproc -use params, only: nanals -use gridio_efsoi - -! Fill in documentation here -implicit none -private -public :: scatter_chunks_ob_impact -real(r_kind),public, allocatable, dimension(:,:,:) :: anal_chunk, anal_chunk_prior -real(r_single),public, allocatable, dimension(:,:) :: ensmean_chunk, ensmean_chunk_prior -real(r_single),public, allocatable, dimension(:,:) :: fcerror_chunk, analmean_chunk -contains - -subroutine scatter_chunks_ob_impact -! distribute chunks from grdin (read in controlvec) according to -! decomposition from load_balance -! AFE use statevec_efsoi, only: grdin,grdin2,grdin3 -use statevec_efsoi, only: grdin,grdin3,grdin5 -use gridinfo_efsoi, only: ncdim -implicit none - -integer(i_kind), allocatable, dimension(:) :: scounts, displs, rcounts -real(r_single), allocatable, dimension(:) :: sendbuf,recvbuf,sendbuf2,recvbuf2 -integer(i_kind) :: np, nn, n, nanal, i, ierr - - -allocate(scounts(0:numproc-1)) -allocate(displs(0:numproc-1)) -allocate(rcounts(0:numproc-1)) -! only IO tasks send any data. -! scounts is number of data elements to send to processor np. -! rcounts is number of data elements to recv from processor np. -! displs is displacement into send array for data to go to proc np -do np=0,numproc-1 - displs(np) = np*npts_max*ncdim -enddo -if (nproc <= nanals-1) then - scounts = npts_max*ncdim -else - scounts = 0 -endif -! displs is also the displacement into recv array for data to go into anal_chunk -! on -! task np. -do np=0,numproc-1 - if (np <= nanals-1) then - rcounts(np) = npts_max*ncdim - else - rcounts(np) = 0 - end if -enddo - -! allocate array to hold pieces of state vector on each proc. -allocate(anal_chunk(nanals,npts_max,ncdim)) -if (nproc == 0) print *,'anal_chunk size = ',size(anal_chunk) - -allocate(anal_chunk_prior(nanals,npts_max,ncdim)) -allocate(ensmean_chunk(npts_max,ncdim)) -allocate(ensmean_chunk_prior(npts_max,ncdim)) -ensmean_chunk = 0. -allocate(sendbuf(numproc*npts_max*ncdim)) -allocate(recvbuf(numproc*npts_max*ncdim)) - -! send and receive buffers. -if (nproc <= nanals-1) then - ! fill up send buffer. - do np=1,numproc - do nn=1,ncdim - do i=1,numptsperproc(np) - n = ((np-1)*ncdim + (nn-1))*npts_max + i - sendbuf(n) = grdin(indxproc(np,i),nn) - enddo - enddo - enddo -end if - -call mpi_alltoallv(sendbuf, scounts, displs, mpi_real4, recvbuf,rcounts, displs,& - mpi_real4, mpi_comm_world, ierr) - -!==> compute ensemble of first guesses on each task, remove mean from anal. -!$omp parallel do schedule(dynamic,1) private(nn,i,nanal,n) -do nn=1,ncdim - do i=1,numptsperproc(nproc+1) - - do nanal=1,nanals - n = ((nanal-1)*ncdim + (nn-1))*npts_max + i - anal_chunk(nanal,i,nn) = recvbuf(n) - enddo - - ensmean_chunk(i,nn) = sum(anal_chunk(:,i,nn))/float(nanals) - ensmean_chunk_prior(i,nn) = ensmean_chunk(i,nn) - - ! remove mean from ensemble. - do nanal=1,nanals - anal_chunk(nanal,i,nn) = anal_chunk(nanal,i,nn)-ensmean_chunk(i,nn) - anal_chunk_prior(nanal,i,nn)=anal_chunk(nanal,i,nn) - end do - - end do -end do -!$omp end parallel do - -deallocate(sendbuf, recvbuf) - -! Get forecast error at the evaluation time and distribute them to all -! processors -allocate(fcerror_chunk(npts_max,ncdim)) -allocate(analmean_chunk(npts_max,ncdim)) -allocate(sendbuf(numproc*npts_max*ncdim)) -allocate(sendbuf2(numproc*npts_max*ncdim)) -allocate(recvbuf(npts_max*ncdim)) -allocate(recvbuf2(npts_max*ncdim)) -if(nproc == 0) then - ! Assign sendbuf recbuf from state variables - do np=1,numproc - do nn=1,ncdim - do i=1,numptsperproc(np) - n = ((np-1)*ncdim + (nn-1))*npts_max + i -! AFE sendbuf(n) = grdin2(indxproc(np,i),nn) - sendbuf(n) = grdin5(indxproc(np,i),nn) - sendbuf2(n) = grdin3(indxproc(np,i),nn) - end do - end do - end do - call mpi_scatter(sendbuf,ncdim*npts_max,mpi_real4,recvbuf, & - & ncdim*npts_max,mpi_real4,0,mpi_comm_world,ierr) - call mpi_scatter(sendbuf2,ncdim*npts_max,mpi_real4,recvbuf2, & - & ncdim*npts_max,mpi_real4,0,mpi_comm_world,ierr) -else - call mpi_scatter(sendbuf,ncdim*npts_max,mpi_real4,recvbuf, & - & ncdim*npts_max,mpi_real4,0,mpi_comm_world,ierr) - call mpi_scatter(sendbuf2,ncdim*npts_max,mpi_real4,recvbuf2, & - & ncdim*npts_max,mpi_real4,0,mpi_comm_world,ierr) -end if - -!$omp parallel do schedule(dynamic,1) private(nn,i,n) - - do nn=1,ncdim - do i=1,numptsperproc(nproc+1) - n = (nn-1)*npts_max + i - fcerror_chunk(i,nn) = recvbuf2(n) - analmean_chunk(i,nn) = recvbuf(n) - end do - end do - deallocate(sendbuf,recvbuf,sendbuf2,recvbuf2) -! AFE if(allocated(grdin2)) deallocate(grdin2) - if(allocated(grdin5)) deallocate(grdin5) - call divide_weight(analmean_chunk) - call divide_weight(ensmean_chunk(:,:)) - -call destroy_weight() - -end subroutine scatter_chunks_ob_impact - -end module scatter_chunks_efsoi diff --git a/util/EFSOI_Utilities/src/statevec_efsoi.f90 b/util/EFSOI_Utilities/src/statevec_efsoi.f90 deleted file mode 100644 index 5b6bfe26a1..0000000000 --- a/util/EFSOI_Utilities/src/statevec_efsoi.f90 +++ /dev/null @@ -1,297 +0,0 @@ -module statevec_efsoi -!$$$ module documentation block -! -! module: statevec_efsoi read ensemble members, ensemble mean forecasts, -! ensemble mean analysis and verification. -! Assign and compute forecast perturbations -! and ensemble mean forecast errors from -! information read -! -! prgmmr: whitaker org: esrl/psd date: 2009-02-23 -! prgmmr: groff org: emc date: 2018-05-14 -! -! abstract: io needed for efsoi calculations -! -! Public Subroutines: -! init_statevec_efsoi: read anavinfo table for defining EFSOI state vector -! read_state_efsoi: read ensemble members, ensemble mean forecasts, -! ensemble mean analyses and verification. Assign -! and compute forecast perturbations and forecast errors -! from information read. -! statevec_efsoi_cleanup: deallocate allocatable arrays. -! -! Public Variables: -! nanals: (integer scalar) number of ensemble members (from module params) -! nlevs: number of analysis vertical levels (from module params). -! -! nc3d: number of 3D control variables -! nc2d: number of 2D control variables -! cvars3d: names of 3D control variables -! cvars2d: names of 2D control variables -! ncdim: total number of 2D fields to update (nc3d*nlevs+nc2d) -! index_pres: an index array with pressure value for given state variable -! -! Modules Used: mpisetup, params, kinds, gridio, gridinfo_efsoi, mpeu_util, constants -! -! program history log: -! 2009-02-23 Initial version (as statevec). -! 2009-11-28 revamped to improve IO speed -! 2015-06-29 add multiple time levels to background -! 2016-05-02 shlyaeva: Modification for reading state vector from table -! 2016-09-07 shlyaeva: moved distribution of ens members to loadbal -! 2016-11-29 shlyaeva: module renamed to controlvec (from statevec); gridinfo_efsoi -! init and cleanup are called from here now -! 2018-05-14 Groff: Adapted from enkf controlvec.f90 to provide -! io functionality necessary for efsoi calculations -! -! attributes: -! language: f95 -! -!$$$ - -use mpisetup -use gridio_efsoi, only: readgriddata_efsoi, get_weight -use gridinfo_efsoi, only: getgridinfo_efsoi, gridinfo_cleanup_efsoi, & - npts, vars3d_supported, vars2d_supported, ncdim -use params, only: nlevs, fgfileprefixes, reducedgrid, & - nanals, nlons, nlats, read_member_forecasts, & - read_verification, read_ensmean_forecast, & - read_analysis_mean, eft, andataname, & - datehr, forecast_impact, gdatehr -use kinds, only: r_kind, i_kind, r_double, r_single -use mpeu_util, only: gettablesize, gettable, getindex -use constants, only: max_varname_length -implicit none - -private - -public :: read_state_efsoi, statevec_cleanup_efsoi, init_statevec_efsoi -real(r_single), public, allocatable, dimension(:,:) :: grdin, grdin2, grdin3, grdin4, grdin5 - -integer(i_kind), public :: nc2d, nc3d -character(len=max_varname_length), allocatable, dimension(:), public :: cvars3d -character(len=max_varname_length), allocatable, dimension(:), public :: cvars2d -integer(i_kind), public, allocatable, dimension(:) :: index_pres -integer(i_kind), public, allocatable, dimension(:) :: clevels - -contains - -subroutine init_statevec_efsoi() -! read table with state vector variables for efsoi -! (code adapted from GSI state_vectors.f90 init_anasv routine -implicit none -character(len=*),parameter:: rcname='anavinfo' -character(len=*),parameter:: tbname='state_vector_efsoi::' -character(len=256),allocatable,dimension(:):: utable -character(len=20) var,source,funcof -integer(i_kind) luin,ii,i,ntot, k, nvars -integer(i_kind) ilev, itracer - -! load file -luin=914 -open(luin,file=rcname,form='formatted') - -! Scan file for desired table first -! and get size of table -call gettablesize(tbname,luin,ntot,nvars) - -! Get contents of table -allocate(utable(nvars)) -call gettable(tbname,luin,ntot,nvars,utable) - -! release file unit -close(luin) - -! Retrieve each token of interest from table and define -! variables participating in control vector - -! Count variables first -nc2d=0; nc3d=0; ncdim=0; -do ii=1,nvars - read(utable(ii),*) var, ilev, itracer, source, funcof - if(ilev==1) then - nc2d=nc2d+1 - ncdim=ncdim+1 - else - nc3d=nc3d+1 - ncdim=ncdim+ilev - endif -enddo - -allocate(cvars3d(nc3d),cvars2d(nc2d),clevels(0:nc3d)) - -! Now load information from table -nc2d=0;nc3d=0 -clevels = 0 -do ii=1,nvars - read(utable(ii),*) var, ilev, itracer, source, funcof - if(ilev==1) then - nc2d=nc2d+1 - cvars2d(nc2d)=trim(adjustl(var)) - else if (ilev==nlevs .or. ilev==nlevs+1) then - nc3d=nc3d+1 - cvars3d(nc3d) = trim(adjustl(var)) - clevels(nc3d) = ilev + clevels(nc3d-1) - else - if (nproc .eq. 0) print *,'Error: only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev - call stop2(503) - endif -enddo - -deallocate(utable) - -allocate(index_pres(ncdim)) -ii=0 -do i=1,nc3d - do k=1,clevels(i)-clevels(i-1) - ii = ii + 1 - index_pres(ii)=k - end do -end do -do i = 1,nc2d - ii = ii + 1 - index_pres(ii) = nlevs+1 -enddo - -! sanity checks -if (ncdim == 0) then - if (nproc == 0) print *, 'Error: there are no variables to update.' - call stop2(501) -endif - -do i = 1, nc2d - if (getindex(vars2d_supported, cvars2d(i))<0) then - if (nproc .eq. 0) then - print *,'Error: 2D variable ', cvars2d(i), ' is not supported in current version.' - print *,'Supported variables: ', vars2d_supported - endif - call stop2(502) - endif -enddo -do i = 1, nc3d - if (getindex(vars3d_supported, cvars3d(i))<0) then - if (nproc .eq. 0) then - print *,'Error: 3D variable ', cvars3d(i), ' is not supported in current version.' - print *,'Supported variables: ', vars3d_supported - endif - call stop2(502) - endif -enddo - -if (nproc == 0) then - print *, '2D control variables: ', cvars2d - print *, '3D control variables: ', cvars3d - print *, 'Control levels: ', clevels - print *, 'nc2d: ', nc2d, ', nc3d: ', nc3d, ', ncdim: ', ncdim -endif - - -call getgridinfo_efsoi(fgfileprefixes(1), reducedgrid, clevels, nc3d) - - - -! Get grid weights for EFSOI -! calculation and evaluation -call get_weight() - - -end subroutine init_statevec_efsoi - - - -! ==================================================================== -subroutine read_state_efsoi() -! read ensemble members on IO tasks -implicit none -real(r_double) :: t1,t2 -integer(i_kind) :: nanal -integer(i_kind) :: ierr - - - -! must at least nanals tasks allocated. -! this does not fully exit, why? AFE -if (numproc < nanals) then - print *,'need at least nanals =',nanals,'MPI tasks,' - print *,'have numproc=',numproc,', exiting ...' - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) -end if - -if (npts < numproc) then - print *,'cannot allocate more than npts =',npts,'MPI tasks, exiting ...' - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) -end if - -! read in whole control vector on i/o procs - keep in memory -! (needed in write_ensemble) -if (nproc <= nanals-1) then - - allocate(grdin(npts,ncdim)) - allocate(grdin2(npts,ncdim)) - allocate(grdin3(npts,ncdim)) - allocate(grdin4(npts,ncdim)) - allocate(grdin5(npts,ncdim)) - - nanal = nproc + 1 - t1 = mpi_wtime() - - - ! Read ensemble member forecasts needed to obtain - ! the forecast perturbations at evaluation forecast time (EFT) - call readgriddata_efsoi(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,read_member_forecasts,nanal=nanal,ft=eft,hr=datehr) - - - if (nproc == 0) then - t2 = mpi_wtime() - print *,'time in readgridata on root',t2-t1,'secs' - end if - - if (nproc==0) then - ! Read ensemble mean forecast from analysis - if(forecast_impact) call readgriddata_efsoi(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin2, & - read_ensmean_forecast,nanal=0,ft=eft,hr=datehr) - - ! Read ensemble mean forecast from first guess - call readgriddata_efsoi(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin3, & - read_ensmean_forecast,nanal=0,ft=eft+6,hr=gdatehr) - - ! Compute One half the sum of ensemble mean forecast quantities - grdin3 = 0.5_r_kind * (grdin3 + grdin2) - - ! Verification at evaluation time - call readgriddata_efsoi(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin4, & - read_verification,infilename=andataname) - - ! [(0.5*(e_f + e_g)) / (nanals - 1)] - grdin3 = (grdin3 - grdin4) / real(nanals-1,r_kind) - ! Normalize for surface pressure ----- (This needs to be corrected) ----- - ! AFE trying making this like EXP-efso_fv3-CWB (which has the line - ! ommented - note reads: Delete the normalization of surface pressure; - ! The calculation of (e_t|0 + e_t|-6) do not need normalization though - ! dimensional analysis. - !grdin3(:,ncdim) = grdin3(:,ncdim) / grdin4(:,3) - ! Read ensemble mean analysis, used in evolving localization - if(forecast_impact) call readgriddata_efsoi(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin5, & - read_analysis_mean,nanal=0,ft=0,hr=datehr) - end if -end if - -end subroutine read_state_efsoi - -subroutine statevec_cleanup_efsoi() -! deallocate module-level allocatable arrays. -if (allocated(cvars3d)) deallocate(cvars3d) -if (allocated(cvars2d)) deallocate(cvars2d) -if (allocated(clevels)) deallocate(clevels) -if (allocated(index_pres)) deallocate(index_pres) -if (nproc <= nanals-1 .and. allocated(grdin)) deallocate(grdin) -if (nproc <= nanals-1 .and. allocated(grdin2)) deallocate(grdin2) -if (nproc <= nanals-1 .and. allocated(grdin3)) deallocate(grdin3) -if (nproc <= nanals-1 .and. allocated(grdin4)) deallocate(grdin4) -if (nproc <= nanals-1 .and. allocated(grdin5)) deallocate(grdin5) -call gridinfo_cleanup_efsoi() -end subroutine statevec_cleanup_efsoi - -end module statevec_efsoi