From 03265576ebe018c79043c51b5c6cf527861b2bce Mon Sep 17 00:00:00 2001 From: Kyle Gerheiser <3209794+kgerheiser@users.noreply.github.com> Date: Fri, 18 Sep 2020 08:28:52 -0400 Subject: [PATCH 1/8] Update NetCDF and HDF5 modules on Hera NCEPLIBS was built on Hera with a different NetCDF and it wasn't updated in the build module. Users could not run the wgrib2 executable with this NetCDF loaded. --- modulefiles/build.hera | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modulefiles/build.hera b/modulefiles/build.hera index cbae95e56..80be679f4 100644 --- a/modulefiles/build.hera +++ b/modulefiles/build.hera @@ -24,6 +24,6 @@ module load landsfcutil/2.4.1 module load wgrib2/2.0.8 module use /scratch1/NCEPDEV/nems/emc.nemspara/soft/modulefiles -module load hdf5_parallel/1.10.6 -module load netcdf_parallel/4.7.4 +module load hdf5_parallel/1.10.6.release +module load netcdf_parallel/4.7.4.release module load esmf/8.0.0_ParallelNetCDF From e3f82aaff71f2e30ac2d3bce84aac43671c5fc1d Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Fri, 18 Sep 2020 09:56:42 -0400 Subject: [PATCH 2/8] Convert all scripts to bash Convert all scripts to bash, which is standard on most machines. See issue #144 for more details. --- driver_scripts/driver_grid.cray.sh | 6 +- driver_scripts/driver_grid.dell.sh | 6 +- driver_scripts/driver_grid.hera.sh | 6 +- driver_scripts/driver_grid.jet.sh | 6 +- driver_scripts/driver_grid.orion.sh | 6 +- reg_tests/chgres_cube/driver.hera.sh | 2 +- reg_tests/chgres_cube/driver.jet.sh | 2 +- reg_tests/chgres_cube/driver.orion.sh | 2 +- reg_tests/grid_gen/c96.uniform.sh | 2 +- reg_tests/grid_gen/driver.hera.sh | 2 +- reg_tests/grid_gen/driver.jet.sh | 2 +- reg_tests/grid_gen/driver.orion.sh | 2 +- reg_tests/grid_gen/gfdl.regional.sh | 2 +- scripts/exemcsfc_global_sfc_prep.sh.ecf | 2 +- ush/chgres_cube.sh | 2 +- ush/emcsfc_ice_blend.sh | 2 +- ush/emcsfc_snow.sh | 15 +-- ush/fv3gfs_driver_grid.sh | 50 ++++---- ush/fv3gfs_filter_topo.sh | 16 +-- ush/fv3gfs_make_grid.sh | 4 +- ush/fv3gfs_make_lake.sh | 6 +- ush/fv3gfs_make_orog.sh | 161 ++++++++++++------------ ush/global_cycle.sh | 14 +-- ush/global_cycle_driver.sh | 4 +- ush/sfc_climo_gen.sh | 13 +- 25 files changed, 163 insertions(+), 172 deletions(-) diff --git a/driver_scripts/driver_grid.cray.sh b/driver_scripts/driver_grid.cray.sh index 28e68c4c7..f2faa4a78 100755 --- a/driver_scripts/driver_grid.cray.sh +++ b/driver_scripts/driver_grid.cray.sh @@ -49,7 +49,7 @@ # 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 8) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TEMP_DIR - and path to the repository # clone - home_dir. # 9) Submit script: "cat $script | bsub". # 10) All files will be placed in "out_dir". @@ -106,12 +106,12 @@ fi #----------------------------------------------------------------------- # Check paths. # home_dir - location of repository. -# TMPDIR - working directory. +# TEMP_DIR - working directory. # out_dir - where files will be placed upon completion. #----------------------------------------------------------------------- export home_dir=$LS_SUBCWD/.. -export TMPDIR=/gpfs/hps3/stmp/$LOGNAME/fv3_grid.$gtype +export TEMP_DIR=/gpfs/hps3/stmp/$LOGNAME/fv3_grid.$gtype export out_dir=/gpfs/hps3/stmp/$LOGNAME/my_grids #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.dell.sh b/driver_scripts/driver_grid.dell.sh index c6da158fc..d2d6d08d0 100755 --- a/driver_scripts/driver_grid.dell.sh +++ b/driver_scripts/driver_grid.dell.sh @@ -51,7 +51,7 @@ # 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 8) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TEMP_DIR - and path to the repository # clone - home_dir. # 9) Submit script: "cat $script | bsub". # 10) All files will be placed in "out_dir". @@ -108,12 +108,12 @@ fi #----------------------------------------------------------------------- # Check paths. # home_dir - location of repository. -# TMPDIR - working directory. +# TEMP_DIR - working directory. # out_dir - where files will be placed upon completion. #----------------------------------------------------------------------- export home_dir=$LS_SUBCWD/.. -export TMPDIR=/gpfs/dell1/stmp/$LOGNAME/fv3_grid.$gtype +export TEMP_DIR=/gpfs/dell1/stmp/$LOGNAME/fv3_grid.$gtype export out_dir=/gpfs/dell1/stmp/$LOGNAME/my_grids #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh index 75c986aba..e00c7de1e 100755 --- a/driver_scripts/driver_grid.hera.sh +++ b/driver_scripts/driver_grid.hera.sh @@ -49,7 +49,7 @@ # 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 8) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TEMP_DIR - and path to the repository # clone - home_dir. # 9) Submit script: "sbatch $script". # 10) All files will be placed in "out_dir". @@ -108,12 +108,12 @@ fi #----------------------------------------------------------------------- # Check paths. # home_dir - location of repository. -# TMPDIR - working directory. +# TEMP_DIR - working directory. # out_dir - where files will be placed upon completion. #----------------------------------------------------------------------- export home_dir=$SLURM_SUBMIT_DIR/.. -export TMPDIR=/scratch2/NCEPDEV/stmp1/$LOGNAME/fv3_grid.$gtype +export TEMP_DIR=/scratch2/NCEPDEV/stmp1/$LOGNAME/fv3_grid.$gtype export out_dir=/scratch2/NCEPDEV/stmp1/$LOGNAME/my_grids #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh index 650b39820..2df3f01a5 100755 --- a/driver_scripts/driver_grid.jet.sh +++ b/driver_scripts/driver_grid.jet.sh @@ -50,7 +50,7 @@ # 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 8) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TEMP_DIR - and path to the repository # clone - home_dir. # 9) Submit script: "sbatch $script". # 10) All files will be placed in "out_dir". @@ -109,12 +109,12 @@ fi #----------------------------------------------------------------------- # Check paths. # home_dir - location of repository. -# TMPDIR - working directory. +# TEMP_DIR - working directory. # out_dir - where files will be placed upon completion. #----------------------------------------------------------------------- export home_dir=$SLURM_SUBMIT_DIR/.. -export TMPDIR=/lfs4/HFIP/emcda/$LOGNAME/stmp/fv3_grid.$gtype +export TEMP_DIR=/lfs4/HFIP/emcda/$LOGNAME/stmp/fv3_grid.$gtype export out_dir=/lfs4/HFIP/emcda/$LOGNAME/stmp/my_grids #----------------------------------------------------------------------- diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh index 208596da3..9eebd32a7 100755 --- a/driver_scripts/driver_grid.orion.sh +++ b/driver_scripts/driver_grid.orion.sh @@ -49,7 +49,7 @@ # 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 8) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TEMP_DIR - and path to the repository # clone - home_dir. # 9) Submit script: "sbatch $script". # 10) All files will be placed in "out_dir". @@ -108,12 +108,12 @@ fi #----------------------------------------------------------------------- # Check paths. # home_dir - location of repository. -# TMPDIR - working directory. +# TEMP_DIR - working directory. # out_dir - where files will be placed upon completion. #----------------------------------------------------------------------- export home_dir=$SLURM_SUBMIT_DIR/.. -export TMPDIR=/work/noaa/stmp/$LOGNAME/fv3_grid.$gtype +export TEMP_DIR=/work/noaa/stmp/$LOGNAME/fv3_grid.$gtype export out_dir=/work/noaa/stmp/$LOGNAME/my_grids #----------------------------------------------------------------------- diff --git a/reg_tests/chgres_cube/driver.hera.sh b/reg_tests/chgres_cube/driver.hera.sh index fa09e7999..c61b72c87 100755 --- a/reg_tests/chgres_cube/driver.hera.sh +++ b/reg_tests/chgres_cube/driver.hera.sh @@ -123,7 +123,7 @@ TEST8=$(sbatch --parsable --ntasks-per-node=6 --nodes=1 -t 0:05:00 -A $PROJECT_C sbatch --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J chgres_summary -o $LOG_FILE -e $LOG_FILE \ --open-mode=append -q $QUEUE -d afterok:$TEST8 << EOF -#!/bin/sh +#!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/chgres_cube/driver.jet.sh b/reg_tests/chgres_cube/driver.jet.sh index b17353920..a9ee38707 100755 --- a/reg_tests/chgres_cube/driver.jet.sh +++ b/reg_tests/chgres_cube/driver.jet.sh @@ -123,7 +123,7 @@ TEST8=$(sbatch --parsable --partition=xjet --nodes=1 --ntasks-per-node=6 -t 0:05 sbatch --partition=xjet --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J chgres_summary -o $LOG_FILE -e $LOG_FILE \ --open-mode=append -q $QUEUE -d afterok:$TEST8 << EOF -#!/bin/sh +#!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/chgres_cube/driver.orion.sh b/reg_tests/chgres_cube/driver.orion.sh index 3026082e1..479eedb96 100755 --- a/reg_tests/chgres_cube/driver.orion.sh +++ b/reg_tests/chgres_cube/driver.orion.sh @@ -126,7 +126,7 @@ TEST8=$(sbatch --parsable --ntasks-per-node=6 --nodes=2 -t 0:10:00 -A $PROJECT_C sbatch --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J chgres_summary -o $LOG_FILE -e $LOG_FILE \ --open-mode=append -q $QUEUE -d afterok:$TEST8 << EOF -#!/bin/sh +#!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/c96.uniform.sh b/reg_tests/grid_gen/c96.uniform.sh index 73da0389d..346c95dbf 100755 --- a/reg_tests/grid_gen/c96.uniform.sh +++ b/reg_tests/grid_gen/c96.uniform.sh @@ -8,7 +8,7 @@ set -x -export TMPDIR=${WORK_DIR}/c96.uniform.work +export TEMP_DIR=${WORK_DIR}/c96.uniform.work export out_dir=${WORK_DIR}/c96.uniform export res=96 diff --git a/reg_tests/grid_gen/driver.hera.sh b/reg_tests/grid_gen/driver.hera.sh index ad19573fd..f2a8536cf 100755 --- a/reg_tests/grid_gen/driver.hera.sh +++ b/reg_tests/grid_gen/driver.hera.sh @@ -71,6 +71,6 @@ TEST2=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_ sbatch --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ --open-mode=append -q $QUEUE -d afterok:$TEST2 << EOF -#!/bin/sh +#!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/driver.jet.sh b/reg_tests/grid_gen/driver.jet.sh index a340582b2..5c1facea3 100755 --- a/reg_tests/grid_gen/driver.jet.sh +++ b/reg_tests/grid_gen/driver.jet.sh @@ -71,6 +71,6 @@ TEST2=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_ sbatch --partition=xjet --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ --open-mode=append -q $QUEUE -d afterok:$TEST2 << EOF -#!/bin/sh +#!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/driver.orion.sh b/reg_tests/grid_gen/driver.orion.sh index 2ae73ff2d..12e56d21b 100755 --- a/reg_tests/grid_gen/driver.orion.sh +++ b/reg_tests/grid_gen/driver.orion.sh @@ -67,6 +67,6 @@ TEST2=$(sbatch --parsable --ntasks-per-node=24 --nodes=1 -t 0:10:00 -A $PROJECT_ sbatch --nodes=1 -t 0:01:00 -A $PROJECT_CODE -J grid_summary -o $LOG_FILE -e $LOG_FILE \ --open-mode=append -q $QUEUE -d afterok:$TEST2 << EOF -#!/bin/sh +#!/bin/bash grep -a '<<<' $LOG_FILE > $SUM_FILE EOF diff --git a/reg_tests/grid_gen/gfdl.regional.sh b/reg_tests/grid_gen/gfdl.regional.sh index 05eb31e6c..00d5c9a82 100755 --- a/reg_tests/grid_gen/gfdl.regional.sh +++ b/reg_tests/grid_gen/gfdl.regional.sh @@ -8,7 +8,7 @@ set -x -export TMPDIR=${WORK_DIR}/gfdl.regional.work +export TEMP_DIR=${WORK_DIR}/gfdl.regional.work export out_dir=${WORK_DIR}/gfdl.regional export res=96 # global resolution in which grid is embedded. diff --git a/scripts/exemcsfc_global_sfc_prep.sh.ecf b/scripts/exemcsfc_global_sfc_prep.sh.ecf index 87593e0dd..04fc0f429 100755 --- a/scripts/exemcsfc_global_sfc_prep.sh.ecf +++ b/scripts/exemcsfc_global_sfc_prep.sh.ecf @@ -1,4 +1,4 @@ -#!/bin/sh +#!/bin/bash #### UNIX Script Documentation Block ################################### # . . diff --git a/ush/chgres_cube.sh b/ush/chgres_cube.sh index fd036e5c8..66261fa89 100755 --- a/ush/chgres_cube.sh +++ b/ush/chgres_cube.sh @@ -8,7 +8,7 @@ # See comments for variable definitions and setup information. #---------------------------------------------------------------------------- -set -x +set -eux #---------------------------------------------------------------------------- # Resolution of target grid. diff --git a/ush/emcsfc_ice_blend.sh b/ush/emcsfc_ice_blend.sh index 35f15d61f..9bceee83e 100755 --- a/ush/emcsfc_ice_blend.sh +++ b/ush/emcsfc_ice_blend.sh @@ -1,4 +1,4 @@ -#!/bin/sh +#!/bin/bash #### UNIX Script Documentation Block ################################### # . . diff --git a/ush/emcsfc_snow.sh b/ush/emcsfc_snow.sh index b738ba49e..7c4c6c429 100755 --- a/ush/emcsfc_snow.sh +++ b/ush/emcsfc_snow.sh @@ -1,4 +1,4 @@ -#!/bin/ksh +#!/bin/bash #### UNIX Script Documentation Block ################################### # . . @@ -56,6 +56,8 @@ fi # are only used in ncep ops when the "prod_util" module is loaded. #----------------------------------------------------------------------- +jlogfile=${jlogfile:-"jlogfile"} + use_prod_util=`echo $UTILROOT` if ((${#use_prod_util} != 0)); then use_prod_util="true" @@ -160,18 +162,19 @@ $WGRIB2 -Sec0 ${IMS_FILE} 2>&1 | grep "grib1 message" status=$? if (( status == 0 )); then # grib 1 file tempdate=$($WGRIB -v $IMS_FILE | head -1) - typeset -L10 IMSDATE10 - IMSDATE10=${tempdate#*D=} + IMSDATE=${tempdate#*D=} else # grib 2 file tempdate=$($WGRIB2 -t $IMS_FILE | head -1) - typeset -L10 IMSDATE10 - IMSDATE10=${tempdate#*d=} + IMSDATE=${tempdate#*d=} fi +IMSDATE10=$(echo $IMSDATE|cut -c1-10) IMSYEAR=$(echo $IMSDATE10 | cut -c1-4) IMSMONTH=$(echo $IMSDATE10 | cut -c5-6) IMSDAY=$(echo $IMSDATE10 | cut -c7-8) IMSHOUR=0 # emc convention is to use 00Z. +pgmout=${pgmout:-OUTPUT} + if test "$use_prod_util" = "true" ; then . prep_step fi @@ -214,8 +217,6 @@ cat > ./fort.41 << ! / ! -pgmout=${pgmout:-OUTPUT} - eval $SNOW2MDLEXEC >> $pgmout 2> errfile rc2=$? diff --git a/ush/fv3gfs_driver_grid.sh b/ush/fv3gfs_driver_grid.sh index 202a8e6a6..c522516eb 100755 --- a/ush/fv3gfs_driver_grid.sh +++ b/ush/fv3gfs_driver_grid.sh @@ -1,4 +1,4 @@ -#!/bin/sh +#!/bin/bash # #----------------------------------------------------------------------- # Driver script to create a cubic-sphere based model grid. @@ -33,7 +33,7 @@ # ./driver_scripts. #----------------------------------------------------------------------- -set -ax +set -eux export machine=${machine:?} @@ -101,19 +101,19 @@ if [ $add_lake = true ]; then fi fi -export TMPDIR=${TMPDIR:?} +export TEMP_DIR=${TEMP_DIR:?} export out_dir=${out_dir:?} export home_dir=${home_dir:-"$PWD/../"} export script_dir=$home_dir/ush -export exec_dir=$home_dir/exec +export exec_dir=${exec_dir:-"$home_dir/exec"} export topo=$home_dir/fix/fix_orog export NCDUMP=${NCDUMP:-ncdump} -rm -fr $TMPDIR -mkdir -p $TMPDIR -cd $TMPDIR ||exit 8 +rm -fr $TEMP_DIR +mkdir -p $TEMP_DIR +cd $TEMP_DIR ||exit 8 #---------------------------------------------------------------------------------- #---------------------------------------------------------------------------------- @@ -142,18 +142,18 @@ if [ $gtype = uniform ] || [ $gtype = stretch ] || [ $gtype = nest ]; then name=C${res}r${rn}n${refine_ratio}_${title} fi - grid_dir=$TMPDIR/$name/grid - orog_dir=$TMPDIR/$name/orog + export grid_dir=$TEMP_DIR/$name/grid + export orog_dir=$TEMP_DIR/$name/orog out_dir=$out_dir/C${res} mkdir -p $out_dir if [ $gtype = nest ]; then filter_dir=$orog_dir # nested grid topography will be filtered online else - filter_dir=$TMPDIR/$name/filter_topo + filter_dir=$TEMP_DIR/$name/filter_topo fi - rm -rf $TMPDIR/$name + rm -rf $TEMP_DIR/$name mkdir -p $grid_dir $orog_dir $filter_dir set +x @@ -178,18 +178,18 @@ if [ $gtype = uniform ] || [ $gtype = stretch ] || [ $gtype = nest ]; then #---------------------------------------------------------------------------------- if [ $machine = WCOSS_C ]; then - touch $TMPDIR/orog.file1 + touch $TEMP_DIR/orog.file1 tile=1 while [ $tile -le $ntiles ]; do - echo "$script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo $TMPDIR " >>$TMPDIR/orog.file1 + echo "$script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo " >>$TEMP_DIR/orog.file1 tile=$(( $tile + 1 )) done - aprun -j 1 -n 4 -N 4 -d 6 -cc depth cfp $TMPDIR/orog.file1 + aprun -j 1 -n 4 -N 4 -d 6 -cc depth cfp $TEMP_DIR/orog.file1 err=$? if [ $err != 0 ]; then exit $err fi - rm $TMPDIR/orog.file1 + rm $TEMP_DIR/orog.file1 else tile=1 while [ $tile -le $ntiles ]; do @@ -198,7 +198,7 @@ if [ $gtype = uniform ] || [ $gtype = stretch ] || [ $gtype = nest ]; then echo "............ Execute fv3gfs_make_orog.sh for tile $tile .................." echo set -x - $script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo $TMPDIR + $script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo err=$? if [ $err != 0 ]; then exit $err @@ -267,10 +267,10 @@ elif [ $gtype = regional_gfdl ] || [ $gtype = regional_esg ]; then halop1=$(( halo + 1 )) tile=7 name=regional - grid_dir=$TMPDIR/${name}/grid - orog_dir=$TMPDIR/${name}/orog - filter_dir=$orog_dir # nested grid topography will be filtered online - rm -rf $TMPDIR/$name + export grid_dir=$TEMP_DIR/${name}/grid + export orog_dir=$TEMP_DIR/${name}/orog + filter_dir=$TEMP_DIR/$name/filter_topo + rm -rf $TEMP_DIR/$name mkdir -p $grid_dir $orog_dir $filter_dir #---------------------------------------------------------------------------------- @@ -364,17 +364,17 @@ elif [ $gtype = regional_gfdl ] || [ $gtype = regional_esg ]; then echo "Begin orography generation at `date`" if [ $machine = WCOSS_C ]; then - echo "$script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo $TMPDIR " >>$TMPDIR/orog.file1 - aprun -j 1 -n 4 -N 4 -d 6 -cc depth cfp $TMPDIR/orog.file1 + echo "$script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo " >>$TEMP_DIR/orog.file1 + aprun -j 1 -n 4 -N 4 -d 6 -cc depth cfp $TEMP_DIR/orog.file1 err=$? - rm $TMPDIR/orog.file1 + rm $TEMP_DIR/orog.file1 else set +x echo echo "............ Execute fv3gfs_make_orog.sh for tile $tile .................." echo set -x - $script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo $TMPDIR + $script_dir/fv3gfs_make_orog.sh $res $tile $grid_dir $orog_dir $script_dir $topo err=$? if [ $err != 0 ]; then exit $err @@ -464,7 +464,7 @@ fi #------------------------------------------------------------------------------------ #------------------------------------------------------------------------------------ -export WORK_DIR=$TMPDIR/sfcfields +export WORK_DIR=$TEMP_DIR/sfcfields export SAVE_DIR=$out_dir/fix_sfc export BASE_DIR=$home_dir export FIX_FV3=$out_dir diff --git a/ush/fv3gfs_filter_topo.sh b/ush/fv3gfs_filter_topo.sh index 4f0b2e8e6..06f537701 100755 --- a/ush/fv3gfs_filter_topo.sh +++ b/ush/fv3gfs_filter_topo.sh @@ -1,5 +1,5 @@ -#!/bin/ksh -set -ax +#!/bin/bash +set -eux #----------------------------------------------------------------------------------------- # @@ -26,10 +26,10 @@ else stretch=1.0 fi -export res=$1 -export griddir=$2 -export orodir=$3 -export outdir=$4 +res=$1 +griddir=$2 +orodir=$3 +outdir=$4 executable=$exec_dir/filter_topo if [ ! -s $executable ]; then @@ -41,8 +41,8 @@ if [ ! -s $executable ]; then exit 1 fi -export mosaic_grid=C${res}_mosaic.nc -export topo_file=oro.C${res} +mosaic_grid=C${res}_mosaic.nc +topo_file=oro.C${res} if [ ! -s $outdir ]; then mkdir -p $outdir ;fi cd $outdir ||exit 8 diff --git a/ush/fv3gfs_make_grid.sh b/ush/fv3gfs_make_grid.sh index 7a7cff105..4510c0500 100755 --- a/ush/fv3gfs_make_grid.sh +++ b/ush/fv3gfs_make_grid.sh @@ -1,6 +1,6 @@ -#!/bin/ksh +#!/bin/bash -set -ax +set -eux #----------------------------------------------------------------------------------------- # diff --git a/ush/fv3gfs_make_lake.sh b/ush/fv3gfs_make_lake.sh index ad0d81f14..69013168f 100755 --- a/ush/fv3gfs_make_lake.sh +++ b/ush/fv3gfs_make_lake.sh @@ -1,10 +1,10 @@ -#!/bin/ksh +#!/bin/bash echo echo "CREATE AND ADD INLAND, LAKE_STATUS, AND LAKE_DEPTH TO THE OROGRAPHY FILES" echo -set -ax +set -eux outdir=$orog_dir indir=$topo @@ -26,7 +26,7 @@ if [ ! -s $exe_inland ]; then exit 1 fi -workdir=$TMPDIR/C${res}/orog/tiles +workdir=$TEMP_DIR/C${res}/orog/tiles if [ ! -s $workdir ]; then mkdir -p $workdir ;fi # Make lake data - diff --git a/ush/fv3gfs_make_orog.sh b/ush/fv3gfs_make_orog.sh index 532c6407f..5b2bae4fa 100755 --- a/ush/fv3gfs_make_orog.sh +++ b/ush/fv3gfs_make_orog.sh @@ -1,56 +1,54 @@ #!/bin/bash -set -ax + +set -eux nargv=$# -export inorogexist=0 +inorogexist=0 -if [ $nargv -eq 6 ]; then # lat-lon grid - export lonb=$1 - export latb=$2 - export outdir=$3 - export script_dir=$4 - export is_latlon=1 - export orogfile="none" - export hist_dir=$5 - export TMPDIR=$6 - export workdir=$TMPDIR/latlon/orog/latlon_${lonb}x${latb} -elif [ $nargv -eq 7 ]; then # cubed-sphere grid - export res=$1 - export lonb=$1 - export latb=$1 - export tile=$2 - export griddir=$3 - export outdir=$4 - export script_dir=$5 - export is_latlon=0 - export orogfile="none" - export hist_dir=$6 - export TMPDIR=$7 - export workdir=$TMPDIR/C${res}/orog/tile$tile +if [ $nargv -eq 5 ]; then # lat-lon grid + lonb=$1 + latb=$2 + outdir=$3 + script_dir=$4 + is_latlon=1 + orogfile="none" + hist_dir=$5 + workdir=$TEMP_DIR/latlon/orog/latlon_${lonb}x${latb} +elif [ $nargv -eq 6 ]; then # cubed-sphere grid + res=$1 + lonb=$1 + latb=$1 + tile=$2 + griddir=$3 + outdir=$4 + script_dir=$5 + is_latlon=0 + orogfile="none" + hist_dir=$6 + workdir=$TEMP_DIR/C${res}/orog/tile$tile elif [ $nargv -eq 8 ]; then # input your own orography files - export res=$1 - export lonb=$1 - export latb=$1 - export tile=$2 - export griddir=$3 - export outdir=$4 - export is_latlon=0 - export inputorog=$5 - export script_dir=$6 - export orogfile=$inputorog:t - export inorogexist=1 - export hist_dir=$7 - export TMPDIR=$8 - export workdir=$TMPDIR/C${res}/orog/tile$tile + res=$1 + lonb=$1 + latb=$1 + tile=$2 + griddir=$3 + outdir=$4 + is_latlon=0 + inputorog=$5 + script_dir=$6 + orogfile=$inputorog:t + inorogexist=1 + hist_dir=$7 + workdir=$TEMP_DIR/C${res}/orog/tile$tile else - echo "number of arguments must be 5 or 6 for cubic sphere grid and 4 for lat-lon grid" - echo "Usage for cubic sphere grid: $0 resolution tile grid_dir out_dir script_dir hist_dir TMPDIR" + echo "Number of arguments must be 6 for cubic sphere grid" + echo "Usage for cubic sphere grid: $0 resolution tile griddir outdir script_dir hist_dir" exit 1 fi -export indir=$hist_dir -export executable=$exec_dir/orog +indir=$hist_dir +executable=$exec_dir/orog if [ ! -s $executable ]; then echo "executable does not exist" exit 1 @@ -61,18 +59,18 @@ if [ ! -s $outdir ]; then mkdir -p $outdir ;fi #jcap is for Gaussian grid #jcap=`expr $latb - 2 ` -export jcap=0 -export NF1=0 -export NF2=0 -export mtnres=1 -export efac=0 -export blat=0 -export NR=0 +jcap=0 +NF1=0 +NF2=0 +mtnres=1 +efac=0 +blat=0 +NR=0 if [ $is_latlon -eq 1 ]; then - export OUTGRID="none" + OUTGRID="none" else - export OUTGRID="C${res}_grid.tile${tile}.nc" + OUTGRID="C${res}_grid.tile${tile}.nc" fi # Make Orograraphy @@ -82,47 +80,42 @@ echo "outdir = $outdir" echo "indir = $indir" cd $workdir -# export MTN_SLM=${indir}/TOP8M_slm.80I1.asc -# cp ${indir}/a_ocean_mask${lonb}x${latb}.txt fort.25 -# cp /home/z1l/GFS_tools/orog/a_ocean_mask${lonb}x${latb}.txt fort.25 -# cp $MTN_SLM fort.14 - cp ${indir}/thirty.second.antarctic.new.bin fort.15 - cp ${indir}/landcover30.fixed . +cp ${indir}/thirty.second.antarctic.new.bin fort.15 +cp ${indir}/landcover30.fixed . # uncomment next line to use the old gtopo30 data. # cp ${indir}/gtopo30_gg.fine.nh fort.235 # use gmted2020 data. - cp ${indir}/gmted2010.30sec.int fort.235 - if [ $inorogexist -eq 1 ]; then - cp $inputorog . - fi +cp ${indir}/gmted2010.30sec.int fort.235 +if [ $inorogexist -eq 1 ]; then + cp $inputorog . +fi - if [ $is_latlon -eq 0 ]; then - cp ${griddir}/$OUTGRID . - fi - cp $executable . +if [ $is_latlon -eq 0 ]; then + cp ${griddir}/$OUTGRID . +fi +cp $executable . - echo $mtnres $lonb $latb $jcap $NR $NF1 $NF2 $efac $blat > INPS - echo $OUTGRID >> INPS - echo $orogfile >> INPS - cat INPS - time $executable < INPS +echo $mtnres $lonb $latb $jcap $NR $NF1 $NF2 $efac $blat > INPS +echo $OUTGRID >> INPS +echo $orogfile >> INPS +cat INPS +time $executable < INPS - if [ $? -ne 0 ]; then - echo "ERROR in running $executable " - exit 1 +if [ $? -ne 0 ]; then + echo "ERROR in running $executable " + exit 1 +else + if [ $is_latlon -eq 1 ]; then + outfile=oro.${lonb}x${latb}.nc else - if [ $is_latlon -eq 1 ]; then - export outfile=oro.${lonb}x${latb}.nc - else - export outfile=oro.C${res}.tile${tile}.nc - fi - - mv ./out.oro.nc $outdir/$outfile - echo "file $outdir/$outfile is created" - echo "Successfully running $executable " - exit 0 + outfile=oro.C${res}.tile${tile}.nc fi + mv ./out.oro.nc $outdir/$outfile + echo "file $outdir/$outfile is created" + echo "Successfully running $executable " + exit 0 +fi exit diff --git a/ush/global_cycle.sh b/ush/global_cycle.sh index 477253681..b44d548a2 100755 --- a/ush/global_cycle.sh +++ b/ush/global_cycle.sh @@ -1,4 +1,4 @@ -#!/bin/ksh +#!/bin/bash ################################################################################ #### UNIX Script Documentation Block # . . @@ -242,9 +242,9 @@ COMIN=${COMIN:-$(pwd)} COMOUT=${COMOUT:-$(pwd)} # Filenames. -XC=${XC} -PREINP=${PREINP} -SUFINP=${SUFINP} +XC=${XC:-" "} +PREINP=${PREINP:-" "} +SUFINP=${SUFINP:-" "} CYCLEXEC=${CYCLEXEC:-$EXECgfs/global_cycle$XC} CDATE=${CDATE:?} @@ -295,10 +295,10 @@ GSI_FILE=${GSI_FILE:-"NULL"} FNTSFA=${FNTSFA:-${COMIN}/${PREINP}sstgrb${SUFINP}} FNACNA=${FNACNA:-${COMIN}/${PREINP}engicegrb${SUFINP}} FNSNOA=${FNSNOA:-${COMIN}/${PREINP}snogrb${SUFINP}} -export INISCRIPT=${INISCRIPT} +export INISCRIPT=${INISCRIPT:-" "} export ERRSCRIPT=${ERRSCRIPT:-'eval [[ $err = 0 ]]'} -export LOGSCRIPT=${LOGSCRIPT} -export ENDSCRIPT=${ENDSCRIPT} +export LOGSCRIPT=${LOGSCRIPT:-" "} +export ENDSCRIPT=${ENDSCRIPT:-" "} # Other variables. export PGMOUT=${PGMOUT:-${pgmout:-'&1'}} export PGMERR=${PGMERR:-${pgmerr:-'&2'}} diff --git a/ush/global_cycle_driver.sh b/ush/global_cycle_driver.sh index 7a8c2a8ce..ce5997cb1 100755 --- a/ush/global_cycle_driver.sh +++ b/ush/global_cycle_driver.sh @@ -1,6 +1,6 @@ -#!/bin/sh +#!/bin/bash -set -ax +set -eux #------------------------------------------------------------------------------------------------- # Update surface fields in FV3 restart files on the cubed-sphere grid diff --git a/ush/sfc_climo_gen.sh b/ush/sfc_climo_gen.sh index c51aa1135..a68e99118 100755 --- a/ush/sfc_climo_gen.sh +++ b/ush/sfc_climo_gen.sh @@ -1,4 +1,4 @@ -#!/bin/ksh +#!/bin/bash #------------------------------------------------------------------------- # Run sfc_climo_gen program to create surface fixed fields, @@ -12,7 +12,7 @@ # # BASE_DIR Location of your repository. # input_sfc_climo_dir Location of raw input surface climo data -# EXEC_DIR Location of program executable +# exec_dir Location of program executable # FIX_DIR Location of 'grid' and 'orog' files # GRIDTYPE Flag to invoke logic for global nests # and regional grids. Valid values are @@ -25,16 +25,13 @@ # WORK_DIR Temporary working directory #------------------------------------------------------------------------- -set -x - -ulimit -s unlimited -ulimit -a +set -eux res=${res:-96} WORK_DIR=${WORK_DIR:-/scratch3/NCEPDEV/stmp1/$LOGNAME/sfc_climo_gen.C${res}} SAVE_DIR=${SAVE_DIR:-$WORK_DIR} BASE_DIR=${BASE_DIR:?} -EXEC_DIR=${EXEC_DIR:-$BASE_DIR/exec} +exec_dir=${exec_dir:-$BASE_DIR/exec} GRIDTYPE=${GRIDTYPE:-NULL} FIX_FV3=${FIX_FV3:-/scratch4/NCEPDEV/global/save/glopara/git/fv3gfs/fix/fix_fv3_gmted2010/C${res}} input_sfc_climo_dir=${input_sfc_climo_dir:?} @@ -80,7 +77,7 @@ vegetation_greenness_method="bilinear" EOF APRUN_SFC=${APRUN_SFC:-"aprun -j 1 -n 6 -N 6"} -$APRUN_SFC $EXEC_DIR/sfc_climo_gen +$APRUN_SFC $exec_dir/sfc_climo_gen rc=$? From 230b35c8e4886ec1d578d2eb08830a509b7d4478 Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Fri, 18 Sep 2020 10:53:42 -0400 Subject: [PATCH 3/8] Add Mark Iredell's vertical coordinate generator Add Mark's program, a 'readme' and some sample scripts to run it. For details see issue #136. --- sorc/CMakeLists.txt | 1 + sorc/vcoord_gen.fd/CMakeLists.txt | 15 ++ sorc/vcoord_gen.fd/driver.f90 | 13 + sorc/vcoord_gen.fd/matrix_utils.f90 | 158 ++++++++++++ sorc/vcoord_gen.fd/vcoord_gen.f90 | 377 ++++++++++++++++++++++++++++ util/vcoord_gen/readme.md | 71 ++++++ util/vcoord_gen/run.cray.sh | 51 ++++ util/vcoord_gen/run.sh | 40 +++ 8 files changed, 726 insertions(+) create mode 100644 sorc/vcoord_gen.fd/CMakeLists.txt create mode 100644 sorc/vcoord_gen.fd/driver.f90 create mode 100644 sorc/vcoord_gen.fd/matrix_utils.f90 create mode 100644 sorc/vcoord_gen.fd/vcoord_gen.f90 create mode 100644 util/vcoord_gen/readme.md create mode 100755 util/vcoord_gen/run.cray.sh create mode 100755 util/vcoord_gen/run.sh diff --git a/sorc/CMakeLists.txt b/sorc/CMakeLists.txt index 2cbca0d16..92838fa68 100644 --- a/sorc/CMakeLists.txt +++ b/sorc/CMakeLists.txt @@ -17,3 +17,4 @@ add_subdirectory(grid_tools.fd) add_subdirectory(chgres_cube.fd) add_subdirectory(orog_mask_tools.fd) add_subdirectory(sfc_climo_gen.fd) +add_subdirectory(vcoord_gen.fd) diff --git a/sorc/vcoord_gen.fd/CMakeLists.txt b/sorc/vcoord_gen.fd/CMakeLists.txt new file mode 100644 index 000000000..7c55675a0 --- /dev/null +++ b/sorc/vcoord_gen.fd/CMakeLists.txt @@ -0,0 +1,15 @@ +set(fortran_src + driver.f90 + matrix_utils.f90 + vcoord_gen.f90) + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8") +endif() + +set(exe_name vcoord_gen) +add_executable(${exe_name} ${fortran_src}) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${exec_dir}) diff --git a/sorc/vcoord_gen.fd/driver.f90 b/sorc/vcoord_gen.fd/driver.f90 new file mode 100644 index 000000000..39e9c0703 --- /dev/null +++ b/sorc/vcoord_gen.fd/driver.f90 @@ -0,0 +1,13 @@ +program driver + implicit none + integer levs,lupp,k + real pbot,psig,ppre,pupp,ptop,dpbot,dpsig,dppre,dpupp,dptop,pmin + real,allocatable:: ak(:),bk(:) + write(0,*) 'Enter levs,lupp,pbot,psig,ppre,pupp,ptop,dpbot,dpsig,dppre,dpupp,dptop' + read *,levs,lupp,pbot,psig,ppre,pupp,ptop,dpbot,dpsig,dppre,dpupp,dptop + allocate(ak(0:levs),bk(0:levs)) + call vcoord_gen(levs,lupp,pbot,psig,ppre,pupp,ptop,dpbot,dpsig,dppre,dpupp,dptop,pmin,ak,bk) + write(0,*) 'pmin=',pmin + print '(2i6)',2,levs + print '(f12.3,f12.8)',(ak(k),bk(k),k=0,levs) +end program diff --git a/sorc/vcoord_gen.fd/matrix_utils.f90 b/sorc/vcoord_gen.fd/matrix_utils.f90 new file mode 100644 index 000000000..003c9bbbf --- /dev/null +++ b/sorc/vcoord_gen.fd/matrix_utils.f90 @@ -0,0 +1,158 @@ +!------------------------------------------------------------------------------- +!$$$ Subprogram documentation block +! +! Subprogram: ludcmp lower and upper triangular decomposition +! Prgmmr: Iredell Org: W/NP23 Date: 2008-08-01 +! +! Abstract: This subprogram decomposes a matrix into a product of +! lower and upper triangular matrices. +! +! Program history log: +! 2008-08-01 Mark Iredell +! +! Usage: call ludcmp(a,n,np,indx,d) +! Input argument list: +! a real(np,np) matrix (will be overwritten) +! n integer order of matrix +! np integer dimension of matrix +! +! Output argument list: +! a real(np,np) LU-decomposed matrix +! (U is upper part of A, including diagonal; +! L is lower part of A, with 1 as diagonal; +! L*U equals original A after permuting) +! indx integer(n) pivot indices +! (original A rows are permuted in order i with indx(i)) +! d real determinant permutation (1 or -1, or 0 if singular) +! (determinant is output diagonal product times d) +! +! Attributes: +! Language: Fortran 90 +! +!$$$ +subroutine ludcmp(a,n,np,indx,d) + implicit none + integer,intent(in):: n,np + real,intent(inout):: a(np,np) + integer,intent(out):: indx(n) + real,intent(out):: d + integer i,j,k,imax + real aamax,sum,dum + real vv(n) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + d=1 + do i=1,n + aamax=0 + do j=1,n + if(abs(a(i,j))>aamax) aamax=abs(a(i,j)) + enddo + if(aamax==0) then + d=0 + return + endif + vv(i)=1/aamax + enddo + do j=1,n + do i=1,j-1 + sum=a(i,j) + do k=1,i-1 + sum=sum-a(i,k)*a(k,j) + enddo + a(i,j)=sum + enddo + aamax=0. + do i=j,n + sum=a(i,j) + do k=1,j-1 + sum=sum-a(i,k)*a(k,j) + enddo + a(i,j)=sum + dum=vv(i)*abs(sum) + if(dum>=aamax) then + imax=i + aamax=dum + endif + enddo + if (j/=imax)then + do k=1,n + dum=a(imax,k) + a(imax,k)=a(j,k) + a(j,k)=dum + enddo + d=-d + vv(imax)=vv(j) + endif + indx(j)=imax + if(a(j,j)==0) then + d=0 + return + endif + if(j/=n)then + dum=1/a(j,j) + do i=j+1,n + a(i,j)=a(i,j)*dum + enddo + endif + enddo +end subroutine +!------------------------------------------------------------------------------- +!$$$ Subprogram documentation block +! +! Subprogram: lubksb lower and upper triangular back substitution +! Prgmmr: Iredell Org: W/NP23 Date: 2008-08-01 +! +! Abstract: This subprogram back substitutes to solve decomposed +! lower and upper triangular matrices as outputted by ludcmp. +! +! Program history log: +! 2008-08-01 Mark Iredell +! +! Usage: call lubksb(a,n,np,indx,b) +! Input argument list: +! a real(np,np) LU-decomposed matrix +! (from ludcmp) +! n integer order of matrix +! np integer dimension of matrix +! indx integer(n) pivot indices +! (from ludcmp) +! b real(n) rhs vector of linear problem (will be overwritten) +! +! Output argument list: +! b real(n) solution of linear problem +! (original A times output B equals original B) +! +! Attributes: +! Language: Fortran 90 +! +!$$$ +subroutine lubksb(a,n,np,indx,b) + implicit none + integer,intent(in):: n,np + real,intent(in):: a(np,np) + integer,intent(in):: indx(n) + real,intent(inout):: b(n) + integer i,j,ii,ll + real sum +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + ii=0 + do i=1,n + ll=indx(i) + sum=b(ll) + b(ll)=b(i) + if (ii/=0)then + do j=ii,i-1 + sum=sum-a(i,j)*b(j) + enddo + elseif(sum/=0) then + ii=i + endif + b(i)=sum + enddo + do i=n,1,-1 + sum=b(i) + do j=i+1,n + sum=sum-a(i,j)*b(j) + enddo + b(i)=sum/a(i,i) + enddo +end subroutine diff --git a/sorc/vcoord_gen.fd/vcoord_gen.f90 b/sorc/vcoord_gen.fd/vcoord_gen.f90 new file mode 100644 index 000000000..0a2aa2471 --- /dev/null +++ b/sorc/vcoord_gen.fd/vcoord_gen.f90 @@ -0,0 +1,377 @@ +!$$$ Subprogram documentation block +! +! Subprogram: vcoord_gen Generates hybrid coordinate interface profiles +! Prgmmr: Iredell Org: W/NP23 Date: 2008-08-01 +! +! Abstract: This subprogram generates hybrid coordinate interface profiles +! from a few given parameters. The hybrid coordinate is intended to start +! out at the bottom in pure sigma and end up at the top in pure pressure, +! with a smooth transition in between. The pressure thickness is close to +! quadratic in pressure, with maximum thicknesses in the middle of the domain. +! The coordinate pressure will have continuous second derivatives in level. +! +! The hybrid coordinate is returned in terms of vectors AK and BK, where +! the interface pressure is defined as A+B*ps where ps is surface pressure +! Thus A=0 in regions of pure sigma and B=0 in regions of pure sigma. +! At the bottom, A(0)=0 and B(0)=1 so that surface pressure is the bottom +! boundary condition, while at the top, A(levs)=ptop and B(levs)=0 so that +! the constant top pressure (which can be zero) is the top boundary condition. +! +! The procedure for the calculation is described in the remarks section below. +! +! Program history log: +! 2008-08-01 Mark Iredell +! +! Usage: call vcoord_gen(levs,lupp,pbot,psig,ppre,pupp,ptop,& +! dpbot,dpsig,dppre,dpupp,dptop,pmin,ak,bk) +! Input argument list: +! levs integer number of levels +! lupp integer number of levels below pupp +! pbot real nominal surface pressure (Pa) +! psig real nominal pressure where coordinate changes +! from pure sigma (Pa) +! ppre real nominal pressure where coordinate changes +! to pure pressure (Pa) +! pupp real nominal pressure where coordinate changes +! to upper atmospheric profile (Pa) +! ptop real pressure at top (Pa) +! dpbot real coordinate thickness at bottom (Pa) +! dpsig real thickness of zone within which coordinate changes +! to pure sigma (Pa) +! dppre real thickness of zone within which coordinate changes +! to pure pressure (Pa) +! dpupp real coordinate thickness at pupp (Pa) +! dptop real coordinate thickness at top (Pa) +! Output argument list: +! pmin real minimum surface pressure (Pa) +! ak real(0:levs) a coordinate values, bottom to top (Pa) +! bk real(0:levs) b coordinate values, bottom to top () +! +! Subprograms called: +! ludcmp lower and upper triangular decomposition +! lubksb lower and upper triangular back substitution +! +! Attributes: +! Language: Fortran 90 +! +! Remarks: +! Example: Create the operational GFS 64-level hybrid coordinate. +! real(8) pmin,ak(0:64),bk(0:64) +! call vcoord_gen(64,64,100000.,99400.,7000.,0.,0.,500.,1200.,18000.,60.,60.,& +! pmin,ak,bk) +! print '(2i6)',2,64 +! print '(f12.3,f12.8)',(ak(k),bk(k),k=0,64) +! end +! +! Graphical description of parameters and zones: +! ptop --- ----- ---------------------- +! ... dptop +! --- zone U (upper atmos) +! ... +! pupp --- ----- ---------------------- +! ... dpupp +! --- ----- +! ... zone P (pure pressure) +! --- +! ... +! ppre --- ----- ---------------------- +! ... +! --- dppre zone T1 (transition 1) +! ... +! --- ----- ---------------------- +! ... +! --- +! ... zone T2 (transition 2) +! --- +! ... +! --- ----- ---------------------- +! ... +! --- dpsig zone T3 (transition 3) +! ... +! psig --- ----- ---------------------- +! ... +! --- ----- zone S (pure sigma) +! ... dpbot +! pbot --- ----- ---------------------- +! +! Detailed procedure description: +! STEP 1. +! The pressure profile is computed with respect to the given reference +! surface pressure pbot. For this surface pressure, the 'sigma' thicknesses +! dsig are assumed to be proportional to a quadratic polynomial in sigma sig +! with zero intercepts sig1 and sig2 somewhere below and above the model +! domain, respectively. That is, +! dsig ~ (sig2-sig)*(sig-sig1)*dk +! Integrating this differential equation gives +! sig = (sig1*exp(c1*k+c2)+sig2)/(exp(c1*k+c2)+1) +! The required boundary conditions sig(0)=1 and sig(levs)=0 +! fix the proportionality and integration constants c1 and c2. +! The two crossing parameters (sig1 and sig2) are determined +! by two input sigma thickness conditions dsig/dk at the bottom and top +! which are respectively given as dpbot/(pbot-pupp) and dpupp/(pbot-pupp). +! The crossing parameters are computed using Newton-Raphson iteration. +! This procedure fixes the pressure profile for surface pressure pbot. +! STEP 2. +! The pressure profile is computed with respect to a minimum surface pressure. +! This minimum surface pressure pmin is yet to be determined. +! Divide the profile into zones: +! zone U (pure pressure) from pupp to ptop +! zone P (pure pressure) from pupp to ppre +! zone T1 (transition 1) from ppre to ppre+dppre +! zone T2 (transition 2) from ppre+dppre to psig-dpsig +! zone T3 (transition 3) from psig-dpsig to psig +! zone S (pure "sigma") from psig to pmin +! (here sigma=p/ps so that d(ln(p))/dk is horizontally uniform) +! The pressure profile in the pure pressure zone P is set from step 1. +! The pressure thicknesses in zone T1 is set to be quadratic in level k. +! The pressure thicknesses in zone T2 is set to be linear in level k. +! The pressure thicknesses in zone T3 is set to be quadratic in level k. +! The pressure profile in the pure sigma zone S is also set from step 1. +! Thus there are nine unknowns: +! the 3 polynomial coefficients in zone T1 +! the 2 polynomial coefficients in zone T2 +! the 3 polynomial coefficients in zone T3 +! and the 1 minimum surface pressure. +! The nine conditions to determine these unknowns are: +! the thickness and its derivative match at zone P and T1 boundary +! the thickness and its derivative match at zone T1 and T2 boundary +! the thickness and its derivative match at zone T2 and T3 boundary +! the thickness and its derivative match at zone T3 and S boundary +! the sum of the thicknesses of zones T1, T2, T3, and S is pmin-ppre +! The unknowns are computed using standard linear decomposition. +! This procedure fixes the pressure profile for surface pressure pmin. +! STEP 3. +! (Step 3 skipped if lupp=levs, in which case pupp=ptop and dpupp=dptop.) +! The pressure in zone U is assumed to be the exponential of a cubic +! polynomial in level k. The function must match the pressure at pupp, +! as well as the thickness and its derivative there, and the pressure +! at ptop+dptop at the second to top level. The latter 3 conditions +! are determined by using standard linear decomposition. +! STEP 4. +! For an arbitrary surface pressure, the pressure profile is an linear +! combination of the pressure profiles for surface pressures pbot and pmin +! p(psfc)=p(pbot)*(psfc-pmin)/(pbot-pmin)+p(pmin)*(pbot-psfc)/(pbot-pmin) +! from which the hybrid coordinate profiles ak and bk are found such that +! p(psfc)=ak+bk*psfc +! +!$$$ +subroutine vcoord_gen(levs,lupp,pbot,psig,ppre,pupp,ptop,& + dpbot,dpsig,dppre,dpupp,dptop,pmin,ak,bk) + implicit none + integer,intent(in):: levs,lupp + real,intent(in):: pbot,psig,ppre,pupp,ptop + real,intent(in):: dpbot,dpsig,dppre,dpupp,dptop + real,intent(out):: pmin,ak(0:levs),bk(0:levs) + integer,parameter:: lo=100,li=10 ! outer and inner N-R iterations + real pdif ! thickness from pbot to pupp + real delb ! delta sig at bot + real delt ! delta sig at top + real sig1 ! crossing parameter 1 + real sig2 ! crossing parameter 2 + real c1 ! proportionality constant + real c2 ! integration constant + real sig ! sig variable + real dsig ! delta sig variable + real delbio0 ! initial guess at delta sig at bot + real deltio0 ! initial guess at delta sig at top + real delbio ! guess at delta sig at bot + real deltio ! guess at delta sig at top + real c1sig1 ! d(c1)/d(sig1) + real c1sig2 ! d(c1)/d(sig2) + real p(2) ! rhs in N-R iteration + real fjac(2,2) ! lhs in N-R iteration + integer indx(2) ! permutations in N-R iteration + real ppred ! pressure at T1-T2 boundary + real spre ! sig at P-T1 boundary + real spred ! sig at T1-T2 boundary + real rkpre ! level at P-T1 boundary + real rkpred ! level at T1-T2 boundary + real pkpre ! dp/dk at P-T1 boundary + real pkkpre ! d2p/dk2 at P-T1 boundary + real psigd ! pressure at T2-T3 boundary + real ssig ! sig at T3-S boundary + real ssigd ! sig at T2-T3 boundary + real rksig ! level at T3-S boundary + real rksigd ! level at T2-T3 boundary + real pksig ! dp/dk at T3-S boundary + real pkksig ! d2p/dk2 at T3-S boundary + real p2sig ! pressure at T3-S boundary for pmin surface pressure + real p2sigd ! pressure at T2-T3 boundary for pmin surface pressure + real p2pred ! pressure at T1-T2 boundary for pmin surface pressure + real x2(9) ! rhs in linear solver + real a2(9,9) ! lhs in linear solver + integer indx2(9) ! permutations in linear solver + real pkupp ! dp/dk at U-P boundary + real pkkupp ! d2p/dk2 at U-P boundary + real x3(3) ! rhs in linear solver + real a3(3,3) ! lhs in linear solver + integer indx3(3) ! permutations in linear solver + real p1 ! pressure variable for pbot surface pressure + real p2 ! pressure variable for pmin surface pressure + real d ! determinant permutation + integer io,ii,k +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! STEP 1. + pdif=pbot-pupp + delb=dpbot/pdif + delt=dpupp/pdif + sig1=1+delb + sig2=-delt + c1=log(-sig2*(1-sig1)/sig1/(sig2-1))/lupp + c2=log((sig2-1)/(1-sig1)) + sig=1 + dsig=(sig2-sig)*(sig-sig1)*c1/(sig1-sig2) + delbio0=-dsig + sig=0 + dsig=(sig2-sig)*(sig-sig1)*c1/(sig1-sig2) + deltio0=-dsig +! Newton-Raphson iterations + do io=1,lo + delbio=delbio0+(delb-delbio0)*io/lo + deltio=deltio0+(delt-deltio0)*io/lo + do ii=1,li + c1sig1=-1/(sig1*(1-sig1)*lupp) + c1sig2=-1/(sig2*(sig2-1)*lupp) + sig=1 + dsig=(sig2-sig)*(sig-sig1)*c1/(sig1-sig2) + p(1)=-delbio-dsig + fjac(1,1)=((-c1*(sig+sig2)+(sig-sig1)*c1sig1*(sig1+sig2)) & + *(sig2-sig)/(sig1+sig2)**2) + fjac(1,2)=((c1*(sig+sig1)+(sig2-sig)*c1sig2*(sig1+sig2)) & + *(sig-sig1)/(sig1+sig2)**2) + sig=0 + dsig=(sig2-sig)*(sig-sig1)*c1/(sig1-sig2) + p(2)=-deltio-dsig + fjac(2,1)=((-c1*(sig+sig2)+(sig-sig1)*c1sig1*(sig1+sig2)) & + *(sig2-sig)/(sig1+sig2)**2) + fjac(2,2)=((c1*(sig+sig1)+(sig2-sig)*c1sig2*(sig1+sig2)) & + *(sig-sig1)/(sig1+sig2)**2) + call ludcmp(fjac,2,2,indx,d) + call lubksb(fjac,2,2,indx,p) + sig1=sig1+p(1) + sig2=sig2+p(2) + c1=log(-sig2*(1-sig1)/sig1/(sig2-1))/lupp + c2=log((sig2-1)/(1-sig1)) + enddo + enddo +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! STEP 2. +! Compute minimum surface pressure + ppred=ppre+dppre + spre=(ppre-pupp)/pdif + spred=(ppred-pupp)/pdif + rkpre=(log((spre-sig2)/(sig1-spre))-c2)/c1 + rkpred=(log((spred-sig2)/(sig1-spred))-c2)/c1 + pkpre=pdif*c1*(sig2-spre)*(spre-sig1)/(sig1-sig2) + pkkpre=pkpre*c1*(sig2+sig1-2*spre)/(sig1-sig2) + psigd=psig-dpsig + ssig=(psig-pupp)/pdif + ssigd=(psigd-pupp)/pdif + rksig=(log((ssig-sig2)/(sig1-ssig))-c2)/c1 + rksigd=(log((ssigd-sig2)/(sig1-ssigd))-c2)/c1 + pksig=pdif*c1*(sig2-ssig)*(ssig-sig1)/(sig1-sig2) + pkksig=pksig*c1*(sig2+sig1-2*ssig)/(sig1-sig2) + x2=0 + a2=0 + x2(1)=pkpre + a2(1,1)=1 + a2(1,2)=rkpre + a2(1,3)=rkpre**2 + x2(2)=pkkpre + a2(2,2)=1 + a2(2,3)=2*rkpre + a2(3,1)=1 + a2(3,2)=rkpred + a2(3,3)=rkpred**2 + a2(3,4)=-1 + a2(3,5)=-rkpred + a2(4,2)=1 + a2(4,3)=2*rkpred + a2(4,5)=-1 + a2(5,4)=-1 + a2(5,5)=-rksigd + a2(5,6)=1 + a2(5,7)=rksigd + a2(5,8)=rksigd**2 + a2(6,5)=-1 + a2(6,7)=1 + a2(6,8)=2*rksigd + a2(7,6)=1 + a2(7,7)=rksig + a2(7,8)=rksig**2 + a2(7,9)=-pksig/pbot + a2(8,7)=1 + a2(8,8)=2*rksig + a2(8,9)=-pkksig/pbot + x2(9)=ppre + a2(9,1)=(rkpre-rkpred) + a2(9,2)=(rkpre**2-rkpred**2)/2 + a2(9,3)=(rkpre**3-rkpred**3)/3 + a2(9,4)=(rkpred-rksigd) + a2(9,5)=(rkpred**2-rksigd**2)/2 + a2(9,6)=(rksigd-rksig) + a2(9,7)=(rksigd**2-rksig**2)/2 + a2(9,8)=(rksigd**3-rksig**3)/3 + a2(9,9)=psig/pbot + call ludcmp(a2,9,9,indx2,d) + call lubksb(a2,9,9,indx2,x2) + pmin=x2(9) + p2sig=psig/pbot*pmin + p2sigd=p2sig & + +x2(6)*(rksigd-rksig) & + +x2(7)*(rksigd**2-rksig**2)/2 & + +x2(8)*(rksigd**3-rksig**3)/3 + p2pred=p2sigd & + +x2(4)*(rkpred-rksigd) & + +x2(5)*(rkpred**2-rksigd**2)/2 +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! STEP 3. + if(lupp.lt.levs) then + pkupp=pdif*c1*(sig2-0)*(0-sig1)/(sig1-sig2) + pkkupp=pkupp*c1*(sig2+sig1-2*0)/(sig1-sig2) + x3=0 + a3=0 + x3(1)=pkupp + a3(1,1)=pupp + x3(2)=pkkupp*pupp-pkupp**2 + a3(2,2)=pupp**2 + x3(3)=log((ptop+dptop)/pupp) + a3(3,1)=levs-1-lupp + a3(3,2)=(levs-1-lupp)**2/2 + a3(3,3)=(levs-1-lupp)**3/3 + call ludcmp(a3,3,3,indx3,d) + call lubksb(a3,3,3,indx3,x3) + endif +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! STEP 4. +! Compute hybrid interface values + ak(0)=0 + bk(0)=1 + do k=1,levs-1 + if(k.ge.lupp) then + p1=pupp*exp(x3(1)*(k-lupp)+x3(2)*(k-lupp)**2/2+x3(3)*(k-lupp)**3/3) + else + p1=(sig1*exp(c1*k+c2)+sig2)/(exp(c1*k+c2)+1)*pdif+pupp + endif + if(k.ge.rkpre) then + p2=p1 + elseif(k.ge.rkpred) then + p2=p2pred+x2(1)*(k-rkpred) & + +x2(2)*(k**2-rkpred**2)/2 & + +x2(3)*(k**3-rkpred**3)/3 + elseif(k.ge.rksigd) then + p2=p2sigd+x2(4)*(k-rksigd) & + +x2(5)*(k**2-rksigd**2)/2 + elseif(k.ge.rksig) then + p2=p2sig+x2(6)*(k-rksig) & + +x2(7)*(k**2-rksig**2)/2 & + +x2(8)*(k**3-rksig**3)/3 + else + p2=p1/pbot*pmin + endif + ak(k)=(p2*pbot-p1*pmin)/(pbot-pmin) + bk(k)=(p1-p2)/(pbot-pmin) + enddo + ak(levs)=ptop + bk(levs)=0 +end subroutine diff --git a/util/vcoord_gen/readme.md b/util/vcoord_gen/readme.md new file mode 100644 index 000000000..e6943a53a --- /dev/null +++ b/util/vcoord_gen/readme.md @@ -0,0 +1,71 @@ + + This program generates hybrid coordinate interface profiles + from a few given parameters. The hybrid coordinate is intended to start + out at the bottom in pure sigma and end up at the top in pure pressure, + with a smooth transition in between. The pressure thickness is close to + quadratic in pressure, with maximum thicknesses in the middle of the domain. + The coordinate pressure will have continuous second derivatives in level. + + The hybrid coordinate is returned in terms of vectors AK and BK, where + the interface pressure is defined as A+B*ps where ps is surface pressure + Thus A=0 in regions of pure sigma and B=0 in regions of pure sigma. + At the bottom, A(0)=0 and B(0)=1 so that surface pressure is the bottom + boundary condition, while at the top, A(levs)=ptop and B(levs)=0 so that + the constant top pressure (which can be zero) is the top boundary condition. + + Input argument list: + levs integer number of levels + lupp integer number of levels below pupp + pbot real nominal surface pressure (Pa) + psig real nominal pressure where coordinate changes + from pure sigma (Pa) + ppre real nominal pressure where coordinate changes + to pure pressure (Pa) + pupp real nominal pressure where coordinate changes + to upper atmospheric profile (Pa) + ptop real pressure at top (Pa) + dpbot real coordinate thickness at bottom (Pa) + dpsig real thickness of zone within which coordinate changes + to pure sigma (Pa) + dppre real thickness of zone within which coordinate changes + to pure pressure (Pa) + dpupp real coordinate thickness at pupp (Pa) + dptop real coordinate thickness at top (Pa) + + Outputs a text file containing the 'ak' and 'bk' values - + bottom to top. The conversion to pressure is: + + pressure = ak + (bk * pbot) + + Graphical description of parameters and zones: + ptop --- ----- ---------------------- + ... dptop + --- zone U (upper atmos) + ... + pupp --- ----- ---------------------- + ... dpupp + --- ----- + ... zone P (pure pressure) + --- + ... + ppre --- ----- ---------------------- + ... + --- dppre zone T1 (transition 1) + ... + --- ----- ---------------------- + ... + --- + ... zone T2 (transition 2) + --- + ... + --- ----- ---------------------- + ... + --- dpsig zone T3 (transition 3) + ... + psig --- ----- ---------------------- + ... + --- ----- zone S (pure sigma) + ... dpbot + pbot --- ----- ---------------------- + + For details on the procedure, see the program prolog. diff --git a/util/vcoord_gen/run.cray.sh b/util/vcoord_gen/run.cray.sh new file mode 100755 index 000000000..b6f9c7e13 --- /dev/null +++ b/util/vcoord_gen/run.cray.sh @@ -0,0 +1,51 @@ +#!/bin/bash + +#BSUB -W 0:01 +#BSUB -o log +#BSUB -e log +#BSUB -J vcoord +#BSUB -q debug +#BSUB -R "rusage[mem=100]" +#BSUB -P GFS-DEV + +#------------------------------------------------------------------------------- +# +# Generate a hybrid coordinate interface profile on WCOSS-Cray. +# +# Build the repository using the ./build_all.sh script before running. +# +# Output 'ak' and 'bk' values are placed in $outfile. +# +#------------------------------------------------------------------------------- + +set -x + +source ../../sorc/machine-setup.sh > /dev/null 2>&1 +source ../../modulefiles/build.$target +module list + +outfile="./global_hyblev.txt" + +levs=128 # integer number of levels +lupp=88 # integer number of levels below pupp +pbot=100000.0 # real nominal surface pressure (Pa) +psig=99500.0 # real nominal pressure where coordinate changes + # from pure sigma (Pa) +ppre=7000.0 # real nominal pressure where coordinate changes + # to pure pressure (Pa) +pupp=7000.0 # real nominal pressure where coordinate changes + # to upper atmospheric profile (Pa) +ptop=0.0 # real pressure at top (Pa) +dpbot=240.0 # real coordinate thickness at bottom (Pa) +dpsig=1200.0 # real thickness of zone within which coordinate changes + # to pure sigma (Pa) +dppre=18000.0 # real thickness of zone within which coordinate changes + # to pure pressure (Pa) +dpupp=550.0 # real coordinate thickness at pupp (Pa) +dptop=1.0 # real coordinate thickness at top (Pa) + +rm -f $outfile + +echo $levs $lupp $pbot $psig $ppre $pupp $ptop $dpbot $dpsig $dppre $dpupp $dptop | $PWD/../../exec/vcoord_gen > $outfile + +exit diff --git a/util/vcoord_gen/run.sh b/util/vcoord_gen/run.sh new file mode 100755 index 000000000..b0d5c3aff --- /dev/null +++ b/util/vcoord_gen/run.sh @@ -0,0 +1,40 @@ +#!/bin/bash + +#------------------------------------------------------------------------------- +# +# Generate a hybrid coordinate interface profile. On WCOSS-Cray, use +# 'run.cray.sh'. +# +# Build the repository using the ./build_all.sh script before running. +# +# Output 'ak' and 'bk' values are placed in $outfile. +# +#------------------------------------------------------------------------------- + +set -x + +outfile="./global_hyblev.txt" + +levs=128 # integer number of levels +lupp=88 # integer number of levels below pupp +pbot=100000.0 # real nominal surface pressure (Pa) +psig=99500.0 # real nominal pressure where coordinate changes + # from pure sigma (Pa) +ppre=7000.0 # real nominal pressure where coordinate changes + # to pure pressure (Pa) +pupp=7000.0 # real nominal pressure where coordinate changes + # to upper atmospheric profile (Pa) +ptop=0.0 # real pressure at top (Pa) +dpbot=240.0 # real coordinate thickness at bottom (Pa) +dpsig=1200.0 # real thickness of zone within which coordinate changes + # to pure sigma (Pa) +dppre=18000.0 # real thickness of zone within which coordinate changes + # to pure pressure (Pa) +dpupp=550.0 # real coordinate thickness at pupp (Pa) +dptop=1.0 # real coordinate thickness at top (Pa) + +rm -f $outfile + +echo $levs $lupp $pbot $psig $ppre $pupp $ptop $dpbot $dpsig $dppre $dpupp $dptop | ../../exec/vcoord_gen > $outfile + +exit From dbcad06d54c339d0d1667360a9ca48753dd7a645 Mon Sep 17 00:00:00 2001 From: gsketefian <31046882+gsketefian@users.noreply.github.com> Date: Fri, 25 Sep 2020 15:33:54 -0600 Subject: [PATCH 4/8] Update make_solo_mosaic to handle longer file names Increase the length of the string "mosaic_name" in make_solo_mosaic.c from 128 to 255 characters. This is needed for the regional workflow where "mosaic_name" includes the absolute path and can exceed 128 characters. In that case, make_solo_mosaic fails without a clear error message. --- sorc/fre-nctools.fd/tools/make_solo_mosaic/make_solo_mosaic.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/fre-nctools.fd/tools/make_solo_mosaic/make_solo_mosaic.c b/sorc/fre-nctools.fd/tools/make_solo_mosaic/make_solo_mosaic.c index 82f7ee78f..f78ec54a5 100644 --- a/sorc/fre-nctools.fd/tools/make_solo_mosaic/make_solo_mosaic.c +++ b/sorc/fre-nctools.fd/tools/make_solo_mosaic/make_solo_mosaic.c @@ -78,7 +78,7 @@ main (int argc, char *argv[]) int contact_tile1_jstart[MAXCONTACT], contact_tile1_jend[MAXCONTACT]; int contact_tile2_istart[MAXCONTACT], contact_tile2_iend[MAXCONTACT]; int contact_tile2_jstart[MAXCONTACT], contact_tile2_jend[MAXCONTACT]; - char mosaic_name[128] = "solo_mosaic"; + char mosaic_name[STRING] = "solo_mosaic"; char grid_descriptor[128] = ""; int c, i, n, m, l, errflg; From 8bef319960e7169c8a948a9ea1990af78df3df86 Mon Sep 17 00:00:00 2001 From: David Wright Date: Tue, 29 Sep 2020 10:54:19 -0400 Subject: [PATCH 5/8] Use FVCOM lake surface conditions New program fvcom_to_FV3 takes lake surface data generated by FVCOM and replaces the corresponding variables in the surface file created by chgres_cube. Input data to fvcom_to_FV3 needs to be in netcdf format and already horizontally interpolated to the model grid. Some QC is done within the code to create a minimum lake ice value of 15%. The code is strongly based upon work done by Eric James (GSL) and Ming Hu (GSL) to get FVCOM surface conditions into HRRRv4. PR #158. --- sorc/CMakeLists.txt | 1 + sorc/fvcom_tools.fd/CMakeLists.txt | 22 + sorc/fvcom_tools.fd/kinds.f90 | 48 + sorc/fvcom_tools.fd/module_ncio.f90 | 2360 +++++++++++++++++++++++ sorc/fvcom_tools.fd/module_nwp.f90 | 280 +++ sorc/fvcom_tools.fd/module_nwp_base.f90 | 126 ++ sorc/fvcom_tools.fd/process_FVCOM.f90 | 217 +++ sorc/fvcom_tools.fd/readme.md | 39 + 8 files changed, 3093 insertions(+) create mode 100644 sorc/fvcom_tools.fd/CMakeLists.txt create mode 100644 sorc/fvcom_tools.fd/kinds.f90 create mode 100644 sorc/fvcom_tools.fd/module_ncio.f90 create mode 100644 sorc/fvcom_tools.fd/module_nwp.f90 create mode 100644 sorc/fvcom_tools.fd/module_nwp_base.f90 create mode 100755 sorc/fvcom_tools.fd/process_FVCOM.f90 create mode 100644 sorc/fvcom_tools.fd/readme.md diff --git a/sorc/CMakeLists.txt b/sorc/CMakeLists.txt index 92838fa68..0bc698c8d 100644 --- a/sorc/CMakeLists.txt +++ b/sorc/CMakeLists.txt @@ -18,3 +18,4 @@ add_subdirectory(chgres_cube.fd) add_subdirectory(orog_mask_tools.fd) add_subdirectory(sfc_climo_gen.fd) add_subdirectory(vcoord_gen.fd) +add_subdirectory(fvcom_tools.fd) diff --git a/sorc/fvcom_tools.fd/CMakeLists.txt b/sorc/fvcom_tools.fd/CMakeLists.txt new file mode 100644 index 000000000..7f0bec898 --- /dev/null +++ b/sorc/fvcom_tools.fd/CMakeLists.txt @@ -0,0 +1,22 @@ +set(fortran_src + kinds.f90 + module_ncio.f90 + module_nwp_base.f90 + module_nwp.f90 + process_FVCOM.f90) + + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -convert big_endian") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU|Clang|AppleClang)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8 -fconvert=big-endian") +endif() + +set(exe_name fvcom_to_FV3) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries( + ${exe_name} + MPI::MPI_Fortran + NetCDF::NetCDF_Fortran) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${exec_dir}) diff --git a/sorc/fvcom_tools.fd/kinds.f90 b/sorc/fvcom_tools.fd/kinds.f90 new file mode 100644 index 000000000..d10a7d797 --- /dev/null +++ b/sorc/fvcom_tools.fd/kinds.f90 @@ -0,0 +1,48 @@ +module kinds +!$$$ module documentation block +! . . . . +! module: kinds +! abstract: Module to hold specification kinds for variable declaration. +! This module is based on (copied from) Paul vanDelst's +! type_kinds module found in the community radiative transfer +! model +! +! module history log: +! +! Subroutines Included: +! +! Functions Included: +! +! remarks: +! The numerical data types defined in this module are: +! r_single - specification kind for single precision (4-byte) real variable +! i_kind - generic specification kind for default integer +! r_kind - generic specification kind for default floating point +! +! +! attributes: +! language: f90 +! +!$$$ end documentation block + implicit none + private +! +! for name string + integer, parameter, public :: len_sta_name = 8 + +! Integer type definitions below + +! Integer types + integer, parameter, public :: i_kind = 4 + integer, parameter, public :: i_short = 2 + integer, parameter, public :: i_byte = 1 +! Real types + integer, parameter, public :: r_single = 4 ! single precision + integer, parameter, public :: r_kind = 8 + +! + real(r_single),parameter,public :: rmissing=-99999.0 + real(i_kind),parameter,public :: imissing=-99999 + real(r_kind),parameter,public :: drmissing=-99999.0 + +end module kinds diff --git a/sorc/fvcom_tools.fd/module_ncio.f90 b/sorc/fvcom_tools.fd/module_ncio.f90 new file mode 100644 index 000000000..5f118ad93 --- /dev/null +++ b/sorc/fvcom_tools.fd/module_ncio.f90 @@ -0,0 +1,2360 @@ +module module_ncio +! +! module: functions to read and write netcdf files +! +! Ming Hu +! +! program history log: +! 2017-11-01 Hu initial build +! +! Subroutines Included: +! + + use netcdf + implicit none + + public :: ncio +! set default to private + private +! + type :: ncio + character(len=256) :: filename + integer :: ncid, status + integer :: debug_level + + integer :: nDims + integer :: ends(4) + integer :: xtype + character(len=40) :: dimname(4) + contains + procedure :: open => open_nc + procedure :: close => close_nc + +! read in dimension from the nc file + procedure :: get_dim => get_dim_nc + +! read in attribute from the nc file + generic :: get_att => get_att_nc_int,get_att_nc_real,get_att_nc_string + procedure :: get_att_nc_int + procedure :: get_att_nc_real + procedure :: get_att_nc_string + +! read in a 1d, 2d, 3d, or 4d field from the nc file + generic :: get_var => get_var_nc_double_1d, get_var_nc_double_2d, & + get_var_nc_double_3d, & + get_var_nc_real_1d,get_var_nc_real_2d, & + get_var_nc_real_3d, & + get_var_nc_short_1d,get_var_nc_short_2d, & + get_var_nc_int_1d,get_var_nc_int_2d, & + get_var_nc_int_3d, & + get_var_nc_char_1d,get_var_nc_char_2d, & + get_var_nc_char_3d + procedure :: get_var_nc_short + procedure :: get_var_nc_short_1d + procedure :: get_var_nc_short_2d + procedure :: get_var_nc_int + procedure :: get_var_nc_int_1d + procedure :: get_var_nc_int_2d + procedure :: get_var_nc_int_3d + procedure :: get_var_nc_real + procedure :: get_var_nc_real_1d + procedure :: get_var_nc_real_2d + procedure :: get_var_nc_real_3d + procedure :: get_var_nc_double + procedure :: get_var_nc_double_1d + procedure :: get_var_nc_double_2d + procedure :: get_var_nc_double_3d + procedure :: get_var_nc_char + procedure :: get_var_nc_char_1d + procedure :: get_var_nc_char_2d + procedure :: get_var_nc_char_3d + +! replace 1d, 2d, 3d, or 4d field from the nc file + generic :: replace_var => replace_var_nc_double_1d, replace_var_nc_double_2d, & + replace_var_nc_double_3d, & + replace_var_nc_real_1d,replace_var_nc_real_2d, & + replace_var_nc_real_3d, & + replace_var_nc_int_1d,replace_var_nc_int_2d, & + replace_var_nc_int_3d, & + replace_var_nc_char_1d,replace_var_nc_char_2d, & + replace_var_nc_char_3d + procedure :: replace_var_nc_int + procedure :: replace_var_nc_int_1d + procedure :: replace_var_nc_int_2d + procedure :: replace_var_nc_int_3d + procedure :: replace_var_nc_real + procedure :: replace_var_nc_real_1d + procedure :: replace_var_nc_real_2d + procedure :: replace_var_nc_real_3d + procedure :: replace_var_nc_double + procedure :: replace_var_nc_double_1d + procedure :: replace_var_nc_double_2d + procedure :: replace_var_nc_double_3d + procedure :: replace_var_nc_char + procedure :: replace_var_nc_char_1d + procedure :: replace_var_nc_char_2d + procedure :: replace_var_nc_char_3d + + procedure :: handle_err + + procedure :: convert_theta2t_2dgrid +!Add a new 3d variable to output file (David.M.Wright) + procedure :: add_new_var => add_new_var_3d + end type ncio + +contains + +subroutine open_nc(this,filename,action,debug_level) +! +! open a netcdf file, set initial debug level +! +! prgmmr: Ming Hu org: GSD date: 2017-11-01 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: filename + character(len=1),intent(in) :: action + integer,intent(in),optional :: debug_level + + integer :: ncid, status + + this%debug_level=20 + if(present(debug_level)) this%debug_level=debug_level + + this%filename=trim(filename) +! open existing netCDF dataset + if(action=="r" .or. action=="R") then + status = nf90_open(path = trim(filename), mode = nf90_nowrite, ncid = ncid) + elseif(action=="w" .or. action=="W") then + status = nf90_open(path = trim(filename), mode = nf90_write, ncid = ncid) + else + write(*,*) 'unknow action :', action + stop 123 + endif + if (status /= nf90_noerr) call this%handle_err(status) + this%ncid=ncid + + if(this%debug_level>0) then + write(*,*) '>>> open file: ',trim(this%filename) + endif + +end subroutine open_nc + +subroutine close_nc(this) +! +! initial netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-04-10 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + + integer :: ncid, status + + ncid=this%ncid +! +! close netCDF dataset + status = nf90_close(ncid) + if (status /= nf90_noerr) call this%handle_err(status) + +end subroutine close_nc + +subroutine get_att_nc_real(this,attname,rval) +! +! get attribute in wrf netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-10-04 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: attname + real, intent(out) :: rval + + integer :: ncid, status + +! open existing netCDF dataset + ncid=this%ncid + +! get date from exisiting NC file + status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), rval) + if (status /= nf90_noerr) call this%handle_err(status) +! +end subroutine get_att_nc_real + +subroutine get_att_nc_int(this,attname,ival) +! +! get attribute in wrf netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-10-04 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: attname + integer, intent(out) :: ival + + integer :: ncid, status + +! open existing netCDF dataset + ncid=this%ncid + +! get date from exisiting NC file + status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), ival) + if (status /= nf90_noerr) call this%handle_err(status) +! +end subroutine get_att_nc_int + +subroutine get_att_nc_string(this,attname,string) +! +! get attribute in wrf netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-10-04 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: attname + character(len=*), intent(out) :: string + + integer :: ncid, status + +! open existing netCDF dataset + ncid=this%ncid + +! get date from exisiting NC file + status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), string) + if (status /= nf90_noerr) call this%handle_err(status) +! +end subroutine get_att_nc_string + + +subroutine get_dim_nc(this,dimname,dimvalue) +! +! get dimensions in netcdf file +! +! prgmmr: Ming Hu org: GSD/AMB date: 2017-11-01 +! +! abstract: +! + + implicit none +! + class(ncio) :: this + character(len=*), intent(in) :: dimname + integer,intent(out) :: dimvalue + + integer :: ncid, status + integer :: DimId + +! open existing netCDF dataset + ncid=this%ncid + +! get dimension from exisiting NC file + status = nf90_inq_dimid(ncid,trim(dimname), DimId) + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_Inquire_Dimension(ncid, DimId, len = dimvalue) + if (status /= nf90_noerr) call this%handle_err(status) +! +end subroutine get_dim_nc + +!==========================begin replace_var ========================== +subroutine replace_var_nc_char_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + character, intent(in) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='replace_var_nc_char_1d' +! + integer :: i +! +! + ilength=nd1 +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + + call this%replace_var_nc_char(varname,ilength,field) +! +end subroutine replace_var_nc_char_1d + +subroutine replace_var_nc_char_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + character, intent(in) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + character,allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_char_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) field(1,1) + endif +! + call this%replace_var_nc_char(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_char_2d + +subroutine replace_var_nc_char_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + character, intent(in) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + character,allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_char_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j,k) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) field(1,1,1) + endif + + call this%replace_var_nc_char(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_char_3d +! +subroutine replace_var_nc_char(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + character, intent(in) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='replace_var_nc_char' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_CHAR) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_put_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>replace variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine replace_var_nc_char +!--- replace_var_nc_char + +!---- replace real +subroutine replace_var_nc_real_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + real(4), intent(in) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='replace_var_nc_real_1d' +! + integer :: i +! +! + ilength=nd1 +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif +! + call this%replace_var_nc_real(varname,ilength,field) +! +end subroutine replace_var_nc_real_1d + +subroutine replace_var_nc_real_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + real(4), intent(in) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + real(4),allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_real_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + + call this%replace_var_nc_real(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_real_2d + +subroutine replace_var_nc_real_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + real(4), intent(in) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + real(4),allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_real_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j,k) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + + call this%replace_var_nc_real(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_real_3d + +subroutine replace_var_nc_real(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + real(4), intent(in) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='replace_var_nc_real' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_FLOAT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_put_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>replace variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine replace_var_nc_real + +!---- repalce double +subroutine replace_var_nc_double_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + real(8), intent(in) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='replace_var_nc_double_1d' +! + integer :: i +! +! + ilength=nd1 +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif +! + call this%replace_var_nc_double(varname,ilength,field) +! +end subroutine replace_var_nc_double_1d + +subroutine replace_var_nc_double_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + real(8), intent(in) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + real(8),allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_double_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + + call this%replace_var_nc_double(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_double_2d + +subroutine replace_var_nc_double_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + real(8), intent(in) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + real(8),allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_double_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j,k) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + + call this%replace_var_nc_double(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_double_3d +! + +subroutine replace_var_nc_double(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + real(8), intent(in) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='replace_var_nc_double' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_DOUBLE) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_put_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>replace variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine replace_var_nc_double +! +!---- replace int +subroutine replace_var_nc_int_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + integer, intent(in) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_int_1d' +! + integer :: i +! +! + ilength=nd1 +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + + call this%replace_var_nc_int(varname,ilength,field) +! +end subroutine replace_var_nc_int_1d + +subroutine replace_var_nc_int_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + integer, intent(in) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + integer,allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_int_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + + call this%replace_var_nc_int(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_int_2d + +subroutine replace_var_nc_int_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + integer, intent(in) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + integer,allocatable :: temp(:) +! + character*40,parameter :: thissubname='replace_var_nc_int_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + temp(istart:iend)=field(:,j,k) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + + call this%replace_var_nc_int(varname,ilength,temp) + + deallocate(temp) +! +end subroutine replace_var_nc_int_3d +! +subroutine replace_var_nc_int(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + integer, intent(in) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='replace_var_nc_int' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_INT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_put_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>replace variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine replace_var_nc_int +! +!==========================end of replace_var ========================== + +!==========================begin get_var ========================== +subroutine get_var_nc_double_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + real(8), intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_double_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_double(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_double_1d + +subroutine get_var_nc_double_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + real(8), intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + real(8),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_double_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_double(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_double_2d + +subroutine get_var_nc_double_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + real(8), intent(out) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + real(8),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_double_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + call this%get_var_nc_double(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + field(:,j,k)=temp(istart:iend) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3) + endif + deallocate(temp) +! +end subroutine get_var_nc_double_3d +! +subroutine get_var_nc_double(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + real(8), intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_double' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_DOUBLE) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_DOUBLE,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(a,I10)') ' data type : ',this%xtype + write(*,'(a,I10)')' dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(a,I5,I10,2x,a)') ' rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_double + +subroutine get_var_nc_real_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + real(4), intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_real_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_real(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_real_1d + +subroutine get_var_nc_real_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + real(4), intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + real(4),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_real_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_real(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_real_2d + +subroutine get_var_nc_real_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + real(4), intent(out) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + real(4),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_real_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + call this%get_var_nc_real(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + field(:,j,k)=temp(istart:iend) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3) + endif + deallocate(temp) +! +end subroutine get_var_nc_real_3d +! +subroutine get_var_nc_real(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + real(4), intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_real' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_FLOAT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_FLOAT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_real + +subroutine get_var_nc_int_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + integer, intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_int_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_int(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_int_1d + +subroutine get_var_nc_int_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + integer, intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + integer,allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_int_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_int(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_int_2d + +subroutine get_var_nc_int_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + integer, intent(out) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + integer,allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_int_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + call this%get_var_nc_int(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + field(:,j,k)=temp(istart:iend) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + do k=1,nd3 + write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k)) + enddo + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3) + endif + deallocate(temp) +! +end subroutine get_var_nc_int_3d +! +subroutine get_var_nc_int(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + integer, intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_int' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_INT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_int +! +subroutine get_var_nc_short_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + integer(2), intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_short_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_short(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_short_1d +! +subroutine get_var_nc_short_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + integer(2), intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + integer(2),allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_short_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_short(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_short_2d +! +subroutine get_var_nc_short(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + integer(2), intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_short' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_SHORT) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_SHORT,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_short + +subroutine get_var_nc_char_1d(this,varname,nd1,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1 ! size of array dval + character, intent(out) :: field(nd1) ! values of the field read in + integer :: ilength +! + character*40,parameter :: thissubname='get_var_nc_char_1d' +! + integer :: i +! +! + ilength=nd1 + call this%get_var_nc_char(varname,ilength,field) +! + if(nd1==this%ends(1)) then + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) (field(i),i=1,min(nd1,10)) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + endif +! +end subroutine get_var_nc_char_1d + +subroutine get_var_nc_char_2d(this,varname,nd1,nd2,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2 ! size of array dval + character, intent(out) :: field(nd1,nd2) ! values of the field read in + integer :: ilength +! + character,allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_char_2d' +! + integer :: i,j,k + integer :: istart,iend +! +! + ilength=nd1*nd2 + allocate(temp(ilength)) + + call this%get_var_nc_char(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2)) then + do j=1,nd2 + istart=(j-1)*nd1+1 + iend=(j-1)*nd1+nd1 + field(:,j)=temp(istart:iend) + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) field(1,1) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2) + endif + deallocate(temp) +! +end subroutine get_var_nc_char_2d + +subroutine get_var_nc_char_3d(this,varname,nd1,nd2,nd3,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: nd1,nd2,nd3 ! size of array dval + character, intent(out) :: field(nd1,nd2,nd3) ! values of the field read in + integer :: ilength +! + character,allocatable :: temp(:) +! + character*40,parameter :: thissubname='get_var_nc_char_3d' +! + integer :: i,j,k + integer :: length2d + integer :: istart,iend +! +! + length2d=nd1*nd2 + ilength=length2d*nd3 + allocate(temp(ilength)) + + call this%get_var_nc_char(varname,ilength,temp) + + if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then + do k=1,nd3 + do j=1,nd2 + istart=(k-1)*length2d+(j-1)*nd1+1 + iend =(k-1)*length2d+(j-1)*nd1+nd1 + field(:,j,k)=temp(istart:iend) + enddo + enddo +! + if(this%debug_level>100) then + write(*,*) trim(thissubname),' show samples:' + write(*,*) field(1,1,1) + endif + else + write(*,*) trim(thissubname),' ERROR: dimension does not match.' + write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3) + endif + deallocate(temp) +! +end subroutine get_var_nc_char_3d +! +subroutine get_var_nc_char(this,varname,ilength,field) +! +! read in one field +! + use netcdf +! + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname ! name of the field to read + integer, intent(in) :: ilength ! size of array dval + character, intent(out) :: field(ilength) ! values of the field read in +! + integer :: ncid +! + integer :: status + integer :: varid + integer :: ends(4),start(4) + + integer :: length4d,length3d,length2d + integer :: nDims,ndim + integer :: dimids(4) + integer :: xtype + character*40 :: dimname + + character*40,parameter :: thissubname='get_var_nc_char' +! + integer :: i,k +! +! + ncid=this%ncid + +! get variable IDs + status = nf90_inq_varid(ncid, trim(varname), VarId) + if(status /= nf90_NoErr) call this%handle_err(status) + +! get dimensions + ends=1 + start=1 + this%ends=1 + + this%dimname=" " +! get variable type + status = nf90_inquire_variable(ncid, VarId, xtype = xtype) + if(status /= nf90_NoErr) call this%handle_err(status) + if(xtype==NF90_CHAR) then + this%xtype=xtype + else + write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_CHAR,' but read in ',xtype + stop 123 + endif + +! get dimension size + status = nf90_inquire_variable(ncid, VarId, ndims = nDims) + if(status /= nf90_NoErr) call this%handle_err(status) + this%ndims=nDims +! + status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims)) + if(status /= nf90_NoErr) call this%handle_err(status) + do i=1,nDims + dimname=" " + status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim) + if (status /= nf90_noerr) call this%handle_err(status) + ends(i)=ndim + this%ends(i)=ends(i) + this%dimname(i)=trim(dimname) + if(this%ends(i) < 1) then + write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i) + stop 1234 + endif + enddo + length2d=ends(1)*ends(2) + length3d=length2d*ends(3) + length4d=length3d*ends(4) + if(ilength .ne. length4d) then + write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d + stop 123 + endif +! + if(nDims <=4 ) then + status = nf90_get_var(ncid, VarId, field, & + start = start(1:4) , & + count = ends(1:4)) + if(status /= nf90_NoErr) call this%handle_err(status) + else + write(*,*) trim(thissubname),'Error: too many dimensions:',nDims + stop 1234 + endif +! + if(this%debug_level>0) then + write(*,'(a,a)') '>>>read in variable: ',trim(varname) + endif + if(this%debug_level>10) then + write(*,'(8x,a,I10)') 'data type : ',this%xtype + write(*,'(8x,a,I10)') 'dimension size: ',this%nDims + do i=1,this%nDims + write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i)) + enddo + endif +! +end subroutine get_var_nc_char +!==========================end of get_var ========================== + +subroutine handle_err(this,status) + use netcdf + implicit none + class(ncio) :: this +! + integer, intent ( in) :: status + if(status /= nf90_noerr) then + print *, trim(nf90_strerror(status)) + stop "Stopped" + end if +end subroutine handle_err + +subroutine convert_theta2t_2dgrid(this,nx,ny,ps,t2) +! convertt theta T to T + implicit none + class(ncio) :: this + + integer :: nx,ny + real, intent(in ) :: ps(nx,ny) + real, intent(inout) :: t2(nx,ny) + + integer :: i,j + real(8) :: rd,cp,rd_over_cp + + + rd = 2.8705e+2_8 + cp = 1.0046e+3_8 ! specific heat of air @pressure (J/kg/K) + rd_over_cp = rd/cp + + do j=1,ny + do i=1,nx + t2(i,j)=t2(i,j)*(ps(i,j)/1000.0)**rd_over_cp - 273.15 + enddo + enddo + +end subroutine convert_theta2t_2dgrid + +subroutine add_new_var_3d(this,varname,dname1,dname2,dname3,lname,units) +! +! prgmmr: David.M.Wright org: UM/GLERL date: 2020-09-01 +! +! abstract: Add a new variable to sfc_data.nc with dimensions +! (Time, yaxis_1, xaxis_1) +! +! Input: varname = Name of variable to be created in netcdf file +! dname1,dname2,dname3 = 1st,2nd, and 3rd dimension names +! lname = long name output for netcdf variable +! units = units to use in netcdf variable +! + use netcdf + implicit none +! + class(ncio) :: this + character(len=*),intent(in) :: varname,dname1,dname2,dname3 & + ,lname,units + integer :: status, ncid, dim1id, dim2id, dim3id, varid + + status = nf90_inq_dimid(this%ncid, dname1, dim1id) + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_inq_dimid(this%ncid, dname2, dim2id) + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_inq_dimid(this%ncid, dname3, dim3id) + if (status /= nf90_noerr) call this%handle_err(status) + + status = nf90_def_var(this%ncid, varname, nf90_double, & + (/ dim1id, dim2id, dim3id /), varid) + if (status /= nf90_noerr) call this%handle_err(status) + + status = nf90_put_att(this%ncid, varid, 'long_name', lname) + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_put_att(this%ncid, varid, 'units', units) + if (status /= nf90_noerr) call this%handle_err(status) + +end subroutine add_new_var_3d + +end module module_ncio diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90 new file mode 100644 index 000000000..40965cdcf --- /dev/null +++ b/sorc/fvcom_tools.fd/module_nwp.f90 @@ -0,0 +1,280 @@ +! module_nwp.f90 +! David Wright +! University of Michigan and GLERL +! 17 Aug 2020 +! +! This module defines FV3LAM and FVCOM forecast data structure and the method to +! read and write observations from and to those data structures. It is used by +! ingest_FVCOM.f90. +! +! This script is strongly based upon Eric James' (ESRL/GSL) work with HRRR/WRF +! to get FVCOM data into the model. + +module module_nwp + + use kinds, only: r_kind, r_single, i_short, rmissing + use module_nwp_base, only: nwpbase +! use module_map_utils, only: map_util + use module_ncio, only: ncio + + implicit none + + public :: fcst_nwp + public :: nwp_type + + private + type :: nwp_type + character(len=6) :: datatype + integer :: numvar, xlat, xlon, xtime + integer :: i_mask, i_sst, i_ice, i_sfcT, i_iceT + character(len=20), allocatable :: varnames(:) + character(len=20), allocatable :: latname + character(len=20), allocatable :: lonname + character(len=20), allocatable :: dimnameEW + character(len=20), allocatable :: dimnameNS + character(len=20), allocatable :: dimnameTIME + + real(r_kind), allocatable :: nwp_mask(:,:,:) + real(r_kind), allocatable :: nwp_sst(:,:,:) + real(r_kind), allocatable :: nwp_ice(:,:,:) + real(r_kind), allocatable :: nwp_sfcT(:,:,:) + real(r_kind), allocatable :: nwp_iceT(:,:,:) + end type nwp_type + + type, extends(nwp_type) :: fcst_nwp + type(nwpbase), pointer :: head => NULL() + type(nwpbase), pointer :: tail => NULL() + contains + procedure :: initial => initial_nwp + procedure :: list_initial => list_initial_nwp + procedure :: read_n => read_nwp + procedure :: finish => finish_nwp + end type fcst_nwp + + type(ncio) :: ncdata +! type(map_util) :: map + + contains + + subroutine initial_nwp(this,itype) + +! This subroutine defines the number of variables and their names for +! each NWP data type. The indices of the variables are +! also defined for later reference. + + class(fcst_nwp) :: this + + character(len=6), intent(in) :: itype + +! FVCOM grid + if (itype==' FVCOM') then + this%datatype = itype + this%numvar = 4 + + this%i_mask = 1 + this%i_sst = 2 + this%i_ice = 3 + this%i_iceT = 4 + this%i_sfcT = 0 + + allocate(this%varnames(this%numvar)) + this%varnames(1) = 'glmask' + this%varnames(2) = 'tsfc' + this%varnames(3) = 'aice' + this%varnames(4) = 'tisfc' + + allocate(this%latname) + allocate(this%lonname) + this%latname = 'lat' + this%lonname = 'lon' + + allocate(this%dimnameEW) + allocate(this%dimnameNS) + allocate(this%dimnameTIME) + this%dimnameEW = 'lon' + this%dimnameNS = 'lat' + this%dimnameTIME = 'Time' + +! FV3LAM grid + + else if (trim(itype)=='FV3LAM') then + this%datatype = itype + this%numvar = 4 + + this%i_mask = 1 + this%i_sst = 2 + this%i_ice = 3 + this%i_iceT = 4 + this%i_sfcT = 0 + + allocate(this%varnames(this%numvar)) + this%varnames(1) = 'slmsk' + this%varnames(2) = 'tsea' + this%varnames(3) = 'fice' + this%varnames(4) = 'tisfc' + + allocate(this%latname) + allocate(this%lonname) + this%latname = 'yaxis_1' + this%lonname = 'xaxis_1' + + allocate(this%dimnameEW) + allocate(this%dimnameNS) + allocate(this%dimnameTIME) + this%dimnameEW = 'xaxis_1' + this%dimnameNS = 'yaxis_1' + this%dimnameTIME = 'Time' + +! If the data type does not match one of the known types, exit. + + else + write(*,*) 'Unknown data type:', itype + stop 1234 + end if + + this%head => NULL() + this%tail => NULL() + + write(*,*) 'Finished initial_nwp' + write(*,*) ' ' + + end subroutine initial_nwp + + subroutine list_initial_nwp(this) + +! This subroutine lists the setup for NWP data that was done by +! the initial_nwp subroutine. + + class(fcst_nwp) :: this + + integer :: k + + write(*,*) 'List initial setup for ', this%datatype + write(*,*) 'number of variables ', this%numvar + write(*,*) 'variable index: mask, sst, ice, sfcT' + write(*,'(15x,10I3)') this%i_mask, this%i_sst, this%i_ice, & + & this%i_sfcT + write(*,*) 'variable name:' + do k=1,this%numvar + write(*,*) k,trim(this%varnames(k)) + enddo + + write(*,*) 'Finished list_initial_nwp' + write(*,*) ' ' + + end subroutine list_initial_nwp + + subroutine read_nwp(this,filename,itype,numlon,numlat,numtimes,time_to_get,mask,sst,ice,sfcT,iceT) + +! This subroutine initializes arrays to receive the NWP data, +! and opens the file and gets the data. + + class(fcst_nwp) :: this + + character(len=5), intent(in) :: itype + character(len=*), intent(in) :: filename + + integer, intent(in) :: time_to_get + integer, intent(inout) :: numlon, numlat, numtimes +! real(r_single), intent(inout) :: mask(:,:), sst(:,:), ice(:,:), sfcT(:,:) + real(r_kind), intent(inout) :: mask(:,:),sst(:,:),ice(:,:),sfcT(:,:) & + ,iceT(:,:) + +! Open the file using module_ncio.f90 code, and find the number of +! lat/lon points + + call ncdata%open(trim(filename),'r',200) + call ncdata%get_dim(this%dimnameEW,this%xlon) + call ncdata%get_dim(this%dimnameNS,this%xlat) + call ncdata%get_dim(this%dimnameTIME,this%xtime) + + write(*,*) 'number of longitudes for file ', filename, this%xlon + numlon = this%xlon + write(*,*) 'number of latitudes for file ', filename, this%xlat + numlat = this%xlat + write(*,*) 'number of times for file ', filename, this%xtime + numtimes = this%xtime + +! Allocate all the arrays to receive data + + allocate(this%nwp_mask(this%xlon,this%xlat,this%xtime)) + allocate(this%nwp_sst(this%xlon,this%xlat,this%xtime)) + allocate(this%nwp_ice(this%xlon,this%xlat,this%xtime)) + allocate(this%nwp_sfcT(this%xlon,this%xlat,this%xtime)) + allocate(this%nwp_iceT(this%xlon,this%xlat,this%xtime)) + +! Get variables from the data file, but only if the variable is +! defined for that data type. + + if (this%i_mask .gt. 0) then + call ncdata%get_var(this%varnames(this%i_mask),this%xlon, & + this%xlat,this%xtime,this%nwp_mask) + mask = this%nwp_mask(:,:,1) + end if + if (this%i_sst .gt. 0) then + call ncdata%get_var(this%varnames(this%i_sst),this%xlon, & + this%xlat,this%xtime,this%nwp_sst) + sst = this%nwp_sst(:,:,time_to_get) + end if + if (this%i_ice .gt. 0) then + call ncdata%get_var(this%varnames(this%i_ice),this%xlon, & + this%xlat,this%xtime,this%nwp_ice) + ice = this%nwp_ice(:,:,time_to_get) + end if + if (this%i_sfcT .gt. 0) then + call ncdata%get_var(this%varnames(this%i_sfcT),this%xlon, & + this%xlat,this%xtime,this%nwp_sfcT) + sfcT = this%nwp_sfcT(:,:,time_to_get) + end if + if (this%i_iceT .gt. 0) then + call ncdata%get_var(this%varnames(this%i_iceT),this%xlon, & + this%xlat,this%xtime,this%nwp_iceT) + iceT = this%nwp_iceT(:,:,time_to_get) + end if + +! Close the netCDF file. + + call ncdata%close + + write(*,*) 'Finished read_nwp' + write(*,*) ' ' + + end subroutine read_nwp + + subroutine finish_nwp(this) + + class(fcst_nwp) :: this + + type(nwpbase), pointer :: thisobs,thisobsnext + + deallocate(this%varnames) + deallocate(this%latname) + deallocate(this%lonname) + deallocate(this%dimnameEW) + deallocate(this%dimnameNS) + deallocate(this%dimnameTIME) + deallocate(this%nwp_mask) + deallocate(this%nwp_sst) + deallocate(this%nwp_ice) + deallocate(this%nwp_sfcT) + deallocate(this%nwp_iceT) + + thisobs => this%head + if(.NOT.associated(thisobs)) then + write(*,*) 'No memory to release' + return + endif + do while(associated(thisobs)) +! write(*,*) 'destroy ==',thisobs%name + + thisobsnext => thisobs%next + call thisobs%destroy() + thisobs => thisobsnext + enddo + + write(*,*) 'Finished finish_nwp' + write(*,*) ' ' + + end subroutine finish_nwp + +end module module_nwp diff --git a/sorc/fvcom_tools.fd/module_nwp_base.f90 b/sorc/fvcom_tools.fd/module_nwp_base.f90 new file mode 100644 index 000000000..d32e841e4 --- /dev/null +++ b/sorc/fvcom_tools.fd/module_nwp_base.f90 @@ -0,0 +1,126 @@ +! module_nwp_base.f90 +! David Wright +! University of Michigan and GLERL +! 17 Aug 2020 +! +! This module defines nwp observation data structure and the method to +! read and write observations from and to those data structures. It is used by +! ingest_FVCOM.f90. +! +! This script is strongly based upon Eric James' (ESRL/GSL) work with HRRR/WRF. +! + +module module_nwp_base + + use kinds, only: r_kind, r_single, rmissing + + implicit none + + public :: nwpbase + public :: nwplocation + + private + +! Define a nwp observation type. + + type nwplocation + real(r_single) :: lon ! stroke longitude + real(r_single) :: lat ! stroke latitiude + end type nwplocation + +! Define a nwp observation type to contain actual data. + + type, extends(nwplocation) :: nwpbase +! HOW DOES THIS POINTER THING WORK? + type(nwpbase), pointer :: next => NULL() + real(r_single) :: time ! observation time + integer :: numvar ! number of variables in this obs type +! real(r_single), allocatable :: obs(:) ! observation value (# numvar) + real(r_kind), allocatable :: obs(:) + logical :: ifquality ! do these obs include quality info? +! GLM has flash_quality_flag + integer, allocatable :: quality(:) ! if so, quality flags + contains + procedure :: list => list_obsbase + procedure :: alloc => alloc_obsbase + procedure :: destroy => destroy_obsbase + end type nwpbase + + contains + + subroutine list_obsbase(this) + +! This subroutine lists the contents of a base nwp observation + + class(nwpbase) :: this + + integer :: i, numvar + +! Write out the lon, lat, and time of the ob + + write(*,'(a,3f10.3)') 'LIGHTNING OB: longitude, latitude, time =', & + this%lon, this%lat, this%time + +! Loop through all variables and print out obs and quality + + numvar = this%numvar + if (numvar >= 1) then +! MULTI-DIMENSIONAL EXAMPLE IN module_obs_base.f90 + write(*,'(a4,10F12.2)') 'obs=', (this%obs(i),i=1,numvar) + if(this%ifquality) & + write(*,'(a4,10I12)') 'qul=', (this%quality(i),i=1,numvar) + else + write(*,*) 'No obs for this location' + endif + + end subroutine list_obsbase + + subroutine alloc_obsbase(this,numvar,ifquality) + +! This subroutine allocates memory for base nwp observation variables +! Input variables: +! numvar : number of variables in this ob type +! itquality: does this observation include quality information? + + class(nwpbase) :: this + + integer, intent(in) :: numvar + + logical, intent(in), optional :: ifquality + + if (numvar >= 1) then + this%numvar = numvar + + if(allocated(this%obs)) deallocate(this%obs) + allocate(this%obs(numvar)) + + this%ifquality=.false. + if(present(ifquality)) this%ifquality = ifquality + if(this%ifquality) allocate(this%quality(numvar)) + + else + write(*,*) 'alloc_obsbase Error: dimension must be larger than 0:', numvar + stop 1234 + endif + + end subroutine alloc_obsbase + + subroutine destroy_obsbase(this) + +! This subroutine releases memory associated with nwp observations + + class(nwpbase) :: this + + this%numvar = 0 + this%time = 0 + + if(allocated(this%obs)) deallocate(this%obs) + + this%ifquality=.false. + if(allocated(this%quality)) deallocate(this%quality) + + this%next => NULL() + + end subroutine destroy_obsbase + +end module module_nwp_base diff --git a/sorc/fvcom_tools.fd/process_FVCOM.f90 b/sorc/fvcom_tools.fd/process_FVCOM.f90 new file mode 100755 index 000000000..da3be6f9b --- /dev/null +++ b/sorc/fvcom_tools.fd/process_FVCOM.f90 @@ -0,0 +1,217 @@ +! process_FVCOM.f90 +! David Wright +! University of Michigan and GLERL +! 17 Aug 2020 +! +! This is the code to put lake surface temperature and aerial ice +! concentration from GLERL-provided FVCOM forecast files (which have +! already been mapped to the FV3-LAM grid) into sfc_data.nc. +! +! This script will take two variables from the command line: +! 1. Name of FV3 sfc data file (e.g. sfc_data.tile7.halo0.nc) +! 2. Name of FVCOM data file (e.g. fvcom.nc) +! +! To run the script, use the following example, modifying file +! names as needed: +! ./fvcom_to_FV3 sfc_data.tile7.halo0.nc fvcom.nc +! +! Code is strongly based upon Eric James' (ESRL/GSL) work +! to update HRRR/WRF Great Lakes' temperature data with FVCOM. +! Code also relies heavily on Ming Hu's ncio module. + +program process_FVCOM + + use mpi + use kinds, only: r_kind, i_kind, r_single + use module_ncio, only: ncio + use module_nwp, only: fcst_nwp + + implicit none + +! MPI variables + integer :: npe, mype, mypeLocal,ierror +! + +! New object-oriented declarations + + type(ncio) :: geo + type(fcst_nwp) :: fcst + +! Grid variables + + character*180 :: geofile + character*2 :: workPath + character*1 :: char1 + + integer :: MAP_PROJ, NLON, NLAT + integer :: fv3lon, fv3lat, fv3times + integer :: lbclon, lbclat, lbctimes + integer :: i, j, t1, t2 + integer :: num_args, ix + + real :: rad2deg = 180.0/3.1415926 + real :: userDX, userDY, CEN_LAT, CEN_LON + real :: userTRUELAT1, userTRUELAT2, MOAD_CEN_LAT, STAND_LON + real :: truelat1, truelat2, stdlon, lat1, lon1, r_earth + real :: knowni, knownj, dx + real :: one, pi, deg2rad + + character(len=180) :: fv3file + character(len=180) :: fvcomfile + character(len=180), dimension(:), allocatable :: args + + real(r_kind), allocatable :: fv3ice(:,:), fv3sst(:,:) + real(r_kind), allocatable :: fv3sfcT(:,:), fv3mask(:,:) + real(r_kind), allocatable :: fv3iceT(:,:) + real(r_kind), allocatable :: lbcice(:,:), lbcsst(:,:) + real(r_kind), allocatable :: lbcsfcT(:,:), lbcmask(:,:) + real(r_kind), allocatable :: lbciceT(:,:) + +! Declare namelists +! SETUP (general control namelist) : + + integer :: update_type + +! namelist/setup/update_type, t2 + +! MPI setup + call MPI_INIT(ierror) + call MPI_COMM_SIZE(mpi_comm_world,npe,ierror) + call MPI_COMM_RANK(mpi_comm_world,mype,ierror) +! +! NCEP LSF has to use all cores allocated to run this application +! but this if check can make sure only one core run through the real code. +! +if(mype==0) then + +! Get file names from command line arguements + num_args = command_argument_count() + allocate(args(num_args)) + + do ix = 1, num_args + call get_command_argument(ix,args(ix)) + end do + + fv3file=trim(args(1)) + write(*,*) trim(fv3file) + fvcomfile=trim(args(2)) + write(*,*) trim(fvcomfile) + + t2 = 1 +! Obtain grid parameters + + workPath='./' + write(geofile,'(a,a)') trim(workPath), trim(fv3file) + write(*,*) 'sfc data file', trim(geofile) + call geo%open(trim(geofile),'r',200) + call geo%get_dim("xaxis_1",NLON) + call geo%get_dim("yaxis_1",NLAT) + write(*,*) 'NLON,NLAT:', NLON, NLAT + call geo%close + + write(*,*) 'Finished reading sfc_data grid information.' + write(*,*) ' ' + +! Allocate variables for I/O + + allocate(fv3ice(nlon,nlat)) + allocate(fv3sfcT(nlon,nlat)) + allocate(fv3sst(nlon,nlat)) + allocate(fv3mask(nlon,nlat)) + allocate(fv3iceT(nlon,nlat)) + + allocate(lbcice(nlon,nlat)) + allocate(lbcsfcT(nlon,nlat)) + allocate(lbcsst(nlon,nlat)) + allocate(lbcmask(nlon,nlat)) + allocate(lbciceT(nlon,nlat)) + +! Read fv3 sfc_data.nc before update + +! fv3file='sfc_data.nc' + t1=1 + + call fcst%initial('FV3LAM') + call fcst%list_initial + call fcst%read_n(trim(fv3file),'FV3LAM',fv3lon,fv3lat,fv3times,t1,fv3mask,fv3sst,fv3ice,fv3sfcT,fv3iceT) + call fcst%finish + + + write(*,*) 'fv3times: ', fv3times + write(*,*) 'time to use: ', t1 + +! Read FVCOM input datasets + +! fvcomfile='fvcom.nc' +! Space infront of ' FVCOM' below is important!! + call fcst%initial(' FVCOM') + call fcst%list_initial + call fcst%read_n(trim(fvcomfile),' FVCOM',lbclon,lbclat,lbctimes,t2,lbcmask,lbcsst,lbcice,lbcsfcT,lbciceT) + call fcst%finish + +! Check that the dimensions match + + if (lbclon .ne. nlon .or. lbclat .ne. nlat) then + write(*,*) 'ERROR: FVCOM/FV3 dimensions do not match:' + write(*,*) 'lbclon: ', lbclon + write(*,*) 'nlon: ', nlon + write(*,*) 'lbclat: ', lbclat + write(*,*) 'nlat: ', nlat + stop 135 + endif + + write(*,*) 'lbctimes: ', lbctimes + write(*,*) 'time to use: ', t2 + +! Update with FVCOM fields and process +! ice cover data. Ice fraction is set +! to a minimum of 15% due to FV3-LAM +! raising any value below 15% to 15%. + + do j=1,nlat + do i=1,nlon + if (lbcmask(i,j) > 0. .and. lbcsst(i,j) .ge. -90.0) then + !If ice fraction below 15%, set to 0 + if (lbcice(i,j) < 0.15) then + lbcice(i,j) = 0.0 + endif + fv3ice(i,j) = lbcice(i,j) + !If ice in FVCOM, but not in FV3-LAM, change to ice + if (lbcice(i,j) > 0. .and. fv3mask(i,j) == 0.) then + fv3mask(i,j) = 2. + endif + !If ice in FV3-LAM and not FVCOM, remove it from FV3-LAM + if (fv3mask(i,j) == 2. .and. lbcice(i,j) == 0.) then + fv3mask(i,j) = 0. + endif + fv3sst(i,j) = lbcsst(i,j) + 273.15 + fv3sfcT(i,j) = lbcsst(i,j) + 273.15 + fv3iceT(i,j) = lbcsst(i,j) + 273.15 + !If ice exists in FVCOM, change ice surface temp + if (lbcice(i,j) > 0.) then + fv3iceT(i,j) = lbciceT(i,j) + 273.15 + end if + endif + enddo + enddo + +! Write out sfc file again + + call geo%open(trim(fv3file),'w',300) + call geo%replace_var("tsea",NLON,NLAT,fv3sst) + call geo%replace_var("fice",NLON,NLAT,fv3ice) + call geo%replace_var("slmsk",NLON,NLAT,fv3mask) + call geo%replace_var("tisfc",NLON,NLAT,fv3iceT) +! Add_New_Var takes names of (Variable,Dim1,Dim2,Dim3,Long_Name,Units) + call geo%add_new_var('glmsk','xaxis_1','yaxis_1','Time','glmsk','none') + call geo%replace_var('glmsk',NLON,NLAT,lbcmask) + call geo%close + + write(6,*) "=== LOWBC UPDATE SUCCESS ===" + +endif ! mype==0 + +call MPI_FINALIZE(ierror) + + +end program process_FVCOM diff --git a/sorc/fvcom_tools.fd/readme.md b/sorc/fvcom_tools.fd/readme.md new file mode 100644 index 000000000..252c7ddf0 --- /dev/null +++ b/sorc/fvcom_tools.fd/readme.md @@ -0,0 +1,39 @@ +**fvcom_to_FV3.exe** + +**Introduction:** + This code replaces lake surface and lake ice temperature along + with aerial ice concentration generated from Great Lakes + Operational Forecast System (GLOFS), an FVCOM-based model, into + sfc_data.nc. + **NOTE** that the variables in the input files must reside on + the same grid. This means data from FVCOM must be horizontally + interpolated to the FV3 grid. This routine will also force a + minimum ice concentration of 15%. If ice concentration is less + than 15% in FVCOM, it will be set to 0% to avoid FV3 from + changing values less than 15% to 15% and generating unrealistic + lake ice temperatures. + +**Library Dependencies:** + Installation depends on the netCDF library and cmake. + +**Running:** + This routine will take two variables from the command line: + 1. Name of FV3 sfc data file (e.g. sfc_data.tile7.halo0.nc) + which is generated from chgres_cube.exe. + 2. Name of FVCOM data file in netcdf format (e.g. fvcom.nc) + + To run the script, use the following example, modifying file + names as needed: + ./fvcom_to_FV3 sfc_data.tile7.halo0.nc fvcom.nc + Output will be to the sfc data file and include lake surface + and lake ice temperature, and lake ice concentration from FVCOM. + + +This routine is *strongly* based upon Eric James' (ESRL/GSL) work + to update HRRR/WRF Great Lakes' temperature data with FVCOM. + It also relies heavily on Ming Hu's (ESRL/GSL) ncio module. + +**For more information, please contact:** + David Wright + University of Michigan and GLERL + dmwright@umich.edu From 745fdbdee648ef8426bd069638db374f0affb1f9 Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Thu, 1 Oct 2020 09:03:52 -0400 Subject: [PATCH 6/8] CHGRES_CUBE GRIB2 Bug Fix Update chgres_cube to create the wgrib2 inventory file on one mpi task, instead of all tasks. When using all tasks, random failures can occur if one task is still writing the file while another is trying to read it. All programs: replace instances of 'include mpif.h' with 'use mpi'. For details, see #157. --- sorc/chgres_cube.fd/atmosphere.F90 | 8 +++---- sorc/chgres_cube.fd/chgres.F90 | 3 +-- sorc/chgres_cube.fd/input_data.F90 | 8 +++---- sorc/chgres_cube.fd/model_grid.F90 | 15 +++++++++---- sorc/chgres_cube.fd/search_util.f90 | 3 +-- sorc/chgres_cube.fd/surface.F90 | 3 +-- sorc/chgres_cube.fd/utils.f90 | 7 +++--- sorc/global_cycle.fd/cycle.f90 | 8 +++---- sorc/global_cycle.fd/read_write_data.f90 | 27 ++++++++++++------------ sorc/sfc_climo_gen.fd/interp.F90 | 6 ++---- sorc/sfc_climo_gen.fd/model_grid.F90 | 3 +-- sorc/sfc_climo_gen.fd/output.f90 | 3 +-- sorc/sfc_climo_gen.fd/program_setup.f90 | 4 ++-- sorc/sfc_climo_gen.fd/search.f90 | 3 +-- sorc/sfc_climo_gen.fd/source_grid.F90 | 3 +-- sorc/sfc_climo_gen.fd/utils.f90 | 7 +++--- 16 files changed, 53 insertions(+), 58 deletions(-) diff --git a/sorc/chgres_cube.fd/atmosphere.F90 b/sorc/chgres_cube.fd/atmosphere.F90 index 03e7e6869..a76ea36b7 100644 --- a/sorc/chgres_cube.fd/atmosphere.F90 +++ b/sorc/chgres_cube.fd/atmosphere.F90 @@ -144,9 +144,9 @@ module atmosphere subroutine atmosphere_driver(localpet) - implicit none + use mpi - include 'mpif.h' + implicit none integer, intent(in) :: localpet @@ -1407,9 +1407,9 @@ SUBROUTINE VINTG ! LANGUAGE: FORTRAN ! ! - IMPLICIT NONE + use mpi - include 'mpif.h' + IMPLICIT NONE REAL(ESMF_KIND_R8), PARAMETER :: DLTDZ=-6.5E-3*287.05/9.80665 REAL(ESMF_KIND_R8), PARAMETER :: DLPVDRT=-2.5E6/461.50 diff --git a/sorc/chgres_cube.fd/chgres.F90 b/sorc/chgres_cube.fd/chgres.F90 index 0d3fcecbd..20d6b94cc 100644 --- a/sorc/chgres_cube.fd/chgres.F90 +++ b/sorc/chgres_cube.fd/chgres.F90 @@ -9,6 +9,7 @@ program chgres ! !------------------------------------------------------------------------- + use mpi use esmf use atmosphere, only : atmosphere_driver @@ -34,8 +35,6 @@ program chgres ! Initialize mpi and esmf environment. !------------------------------------------------------------------------- - include 'mpif.h' - call mpi_init(ierr) print*,"- INITIALIZE ESMF" diff --git a/sorc/chgres_cube.fd/input_data.F90 b/sorc/chgres_cube.fd/input_data.F90 index 6f324d1c3..f5847fe6d 100644 --- a/sorc/chgres_cube.fd/input_data.F90 +++ b/sorc/chgres_cube.fd/input_data.F90 @@ -1728,9 +1728,9 @@ end subroutine read_input_atm_restart_file subroutine read_input_atm_gaussian_netcdf_file(localpet) - implicit none + use mpi - include 'mpif.h' + implicit none integer, intent(in) :: localpet @@ -2111,9 +2111,9 @@ end subroutine read_input_atm_gaussian_netcdf_file subroutine read_input_atm_tiled_history_file(localpet) - implicit none + use mpi - include 'mpif.h' + implicit none integer, intent(in) :: localpet diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 9beb96c3b..84e84a7a3 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -602,6 +602,8 @@ end subroutine define_input_grid_mosaic subroutine define_input_grid_gfs_grib2(localpet, npets) + use mpi + use wgrib2api use program_setup, only : data_dir_input_grid, & @@ -613,7 +615,7 @@ subroutine define_input_grid_gfs_grib2(localpet, npets) character(len=250) :: the_file - integer :: i, j, rc, clb(2), cub(2) + integer :: i, j, rc, clb(2), cub(2), ierr real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) @@ -631,10 +633,15 @@ subroutine define_input_grid_gfs_grib2(localpet, npets) num_tiles_input_grid = 1 the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid - print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) - rc=grb2_mk_inv(the_file,inv_file) - if (rc /=0) call error_handler("OPENING GRIB2 FILE",rc) + if(localpet == 0) then + print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) + rc=grb2_mk_inv(the_file,inv_file) + if (rc /=0) call error_handler("OPENING GRIB2 FILE",rc) + endif +! Wait for localpet 0 to create inventory. + call mpi_barrier(mpi_comm_world, ierr) + rc = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, & lat=lat4, lon=lon4) if (rc /= 1) call error_handler("READING GRIB2 FILE", rc) diff --git a/sorc/chgres_cube.fd/search_util.f90 b/sorc/chgres_cube.fd/search_util.f90 index 0810b0b9b..58bb3cf37 100644 --- a/sorc/chgres_cube.fd/search_util.f90 +++ b/sorc/chgres_cube.fd/search_util.f90 @@ -34,12 +34,11 @@ subroutine search (field, mask, idim, jdim, tile, field_num, latitude) ! future upgrade. !----------------------------------------------------------------------- + use mpi use esmf implicit none - include 'mpif.h' - integer, intent(in) :: idim, jdim, tile, field_num integer(esmf_kind_i8), intent(in) :: mask(idim,jdim) diff --git a/sorc/chgres_cube.fd/surface.F90 b/sorc/chgres_cube.fd/surface.F90 index 6baa2f9bf..ce1b64198 100644 --- a/sorc/chgres_cube.fd/surface.F90 +++ b/sorc/chgres_cube.fd/surface.F90 @@ -242,6 +242,7 @@ end subroutine surface_driver subroutine interp(localpet) + use mpi use esmf use input_data, only : canopy_mc_input_grid, & @@ -301,8 +302,6 @@ subroutine interp(localpet) implicit none - include 'mpif.h' - integer, intent(in) :: localpet integer :: l(1), u(1) diff --git a/sorc/chgres_cube.fd/utils.f90 b/sorc/chgres_cube.fd/utils.f90 index 3a547eede..337083fd3 100644 --- a/sorc/chgres_cube.fd/utils.f90 +++ b/sorc/chgres_cube.fd/utils.f90 @@ -1,8 +1,8 @@ subroutine error_handler(string, rc) - implicit none + use mpi - include 'mpif.h' + implicit none character(len=*), intent(in) :: string @@ -18,6 +18,7 @@ end subroutine error_handler subroutine netcdf_err( err, string ) + use mpi use netcdf implicit none @@ -26,8 +27,6 @@ subroutine netcdf_err( err, string ) character(len=256) :: errmsg integer :: iret - include "mpif.h" - if( err.EQ.NF90_NOERR )return errmsg = NF90_STRERROR(err) print*,'' diff --git a/sorc/global_cycle.fd/cycle.f90 b/sorc/global_cycle.fd/cycle.f90 index 290ab531f..034183336 100644 --- a/sorc/global_cycle.fd/cycle.f90 +++ b/sorc/global_cycle.fd/cycle.f90 @@ -88,10 +88,10 @@ PROGRAM SFC_DRV ! Added mpi directives. !---------------------------------------------------------------------- + use mpi + IMPLICIT NONE ! - include 'mpif.h' - CHARACTER(LEN=3) :: DONST INTEGER :: IDIM, JDIM, LSOIL, LUGB, IY, IM, ID, IH, IALB INTEGER :: ISOT, IVEGSRC, LENSFC, ZSEA1_MM, ZSEA2_MM, IERR @@ -542,9 +542,9 @@ SUBROUTINE ADJUST_NSST(RLA,RLO,SLMSK_TILE,SLMSK_FG_TILE,SKINT_TILE,& SLMSK_GAUS, DTREF_GAUS, & NSST_DATA - IMPLICIT NONE + USE MPI - include 'mpif.h' + IMPLICIT NONE INTEGER, INTENT(IN) :: LENSFC, LSOIL, IDIM, JDIM, MON, DAY diff --git a/sorc/global_cycle.fd/read_write_data.f90 b/sorc/global_cycle.fd/read_write_data.f90 index 371a42668..5ae7bbeae 100644 --- a/sorc/global_cycle.fd/read_write_data.f90 +++ b/sorc/global_cycle.fd/read_write_data.f90 @@ -63,6 +63,8 @@ subroutine write_data(slifcs,tsffcs,snofcs,tg3fcs,zorfcs, & ! let the model compute it. !------------------------------------------------------------------ + use mpi + implicit none integer, intent(in) :: idim, jdim, lensfc, lsoil @@ -118,8 +120,6 @@ subroutine write_data(slifcs,tsffcs,snofcs,tg3fcs,zorfcs, & real(kind=4), allocatable :: lsoil_data(:), x_data(:), y_data(:) real(kind=8), allocatable :: dum2d(:,:), dum3d(:,:,:) - include "mpif.h" - call mpi_comm_rank(mpi_comm_world, myrank, error) write(rankch, '(i3.3)') (myrank+1) @@ -853,9 +853,9 @@ SUBROUTINE READ_LAT_LON_OROG(RLA,RLO,OROG,OROG_UF,& ! THE "GRID" FILE. !-------------------------------------------------------------- - IMPLICIT NONE + USE MPI - include "mpif.h" + IMPLICIT NONE INTEGER, INTENT(IN) :: IDIM, JDIM, IJDIM @@ -985,9 +985,9 @@ SUBROUTINE NETCDF_ERR( ERR, STRING ) ! AND STOP PROCESSING. !-------------------------------------------------------------- - IMPLICIT NONE + USE MPI - include 'mpif.h' + IMPLICIT NONE INTEGER, INTENT(IN) :: ERR CHARACTER(LEN=*), INTENT(IN) :: STRING @@ -1089,9 +1089,9 @@ SUBROUTINE READ_DATA(TSFFCS,SMCFCS,SNOFCS,STCFCS, & ! SELECTED) FOR A SINGLE CUBED-SPHERE TILE. !----------------------------------------------------------------- - IMPLICIT NONE + USE MPI - include "mpif.h" + IMPLICIT NONE INTEGER, INTENT(IN) :: LSOIL, LENSFC @@ -1541,9 +1541,9 @@ subroutine read_tf_clim_grb(file_sst,sst,rlats_sst,rlons_sst,mlat_sst,mlon_sst,m ! language: f90 ! !$$$ - implicit none + use mpi - include "mpif.h" + implicit none ! declare passed variables and arrays character(*) , intent(in ) :: file_sst @@ -1691,9 +1691,9 @@ subroutine get_tf_clm_dim(file_sst,mlat_sst,mlon_sst) ! machine: ibm rs/6000 sp ! !$$$ - implicit none + use mpi - include "mpif.h" + implicit none ! declare passed variables and arrays character(*) , intent(in ) :: file_sst @@ -1844,10 +1844,9 @@ end subroutine get_dim_nc subroutine nc_check(status) + use mpi use netcdf - include "mpif.h" - integer, intent ( in) :: status integer :: ierr diff --git a/sorc/sfc_climo_gen.fd/interp.F90 b/sorc/sfc_climo_gen.fd/interp.F90 index 1e6141708..41b59f916 100644 --- a/sorc/sfc_climo_gen.fd/interp.F90 +++ b/sorc/sfc_climo_gen.fd/interp.F90 @@ -23,11 +23,10 @@ subroutine interp(localpet, method, input_file) use model_grid use source_grid use utils + use mpi implicit none - include 'mpif.h' - character(len=*), intent(in) :: input_file integer :: rc, localpet @@ -294,11 +293,10 @@ subroutine adjust_for_landice(field, vegt, idim, jdim, field_ch) !----------------------------------------------------------------------- use esmf + use mpi implicit none - include 'mpif.h' - character(len=*), intent(in) :: field_ch integer, intent(in) :: idim, jdim diff --git a/sorc/sfc_climo_gen.fd/model_grid.F90 b/sorc/sfc_climo_gen.fd/model_grid.F90 index 056a57130..d86241556 100644 --- a/sorc/sfc_climo_gen.fd/model_grid.F90 +++ b/sorc/sfc_climo_gen.fd/model_grid.F90 @@ -85,11 +85,10 @@ subroutine define_model_grid(localpet, npets) use netcdf use program_setup use utils + use mpi implicit none - include 'mpif.h' - integer, intent(in) :: localpet, npets character(len=500) :: the_file diff --git a/sorc/sfc_climo_gen.fd/output.f90 b/sorc/sfc_climo_gen.fd/output.f90 index 87016e4af..6055404d5 100644 --- a/sorc/sfc_climo_gen.fd/output.f90 +++ b/sorc/sfc_climo_gen.fd/output.f90 @@ -26,6 +26,7 @@ subroutine output(data_one_tile, lat_one_tile, lon_one_tile, i_mdl, j_mdl, & ! time Time period to be output. !-------------------------------------------------------------------------- + use mpi use esmf use netcdf use utils @@ -36,8 +37,6 @@ subroutine output(data_one_tile, lat_one_tile, lon_one_tile, i_mdl, j_mdl, & implicit none - include 'mpif.h' - integer, intent(in) :: i_mdl, j_mdl, tile integer, intent(in) :: record, time, field_idx diff --git a/sorc/sfc_climo_gen.fd/program_setup.f90 b/sorc/sfc_climo_gen.fd/program_setup.f90 index 05e466485..1cb59e6eb 100644 --- a/sorc/sfc_climo_gen.fd/program_setup.f90 +++ b/sorc/sfc_climo_gen.fd/program_setup.f90 @@ -102,9 +102,9 @@ subroutine read_setup_namelist(localpet) ! localpet mpi task number !----------------------------------------------------------------------- - implicit none + use mpi - include 'mpif.h' + implicit none integer, intent(in) :: localpet diff --git a/sorc/sfc_climo_gen.fd/search.f90 b/sorc/sfc_climo_gen.fd/search.f90 index e7446b096..f62cecc17 100644 --- a/sorc/sfc_climo_gen.fd/search.f90 +++ b/sorc/sfc_climo_gen.fd/search.f90 @@ -28,12 +28,11 @@ subroutine search (field, mask, idim, jdim, tile, field_name) ! field field after missing values are replaced !----------------------------------------------------------------------- + use mpi use esmf implicit none - include 'mpif.h' - character(len=*) :: field_name integer, intent(in) :: idim, jdim, tile diff --git a/sorc/sfc_climo_gen.fd/source_grid.F90 b/sorc/sfc_climo_gen.fd/source_grid.F90 index d73e850a1..ccca00b10 100644 --- a/sorc/sfc_climo_gen.fd/source_grid.F90 +++ b/sorc/sfc_climo_gen.fd/source_grid.F90 @@ -79,12 +79,11 @@ subroutine define_source_grid(localpet, npets, input_file) ! !----------------------------------------------------------------------- + use mpi use netcdf implicit none - include 'mpif.h' - character(len=*), intent(in) :: input_file integer, intent(in) :: localpet, npets diff --git a/sorc/sfc_climo_gen.fd/utils.f90 b/sorc/sfc_climo_gen.fd/utils.f90 index 0941737d1..fe18a17dc 100644 --- a/sorc/sfc_climo_gen.fd/utils.f90 +++ b/sorc/sfc_climo_gen.fd/utils.f90 @@ -9,6 +9,7 @@ module utils subroutine netcdf_err( err, string ) + use mpi use netcdf implicit none @@ -17,8 +18,6 @@ subroutine netcdf_err( err, string ) character(len=256) :: errmsg integer :: ierr - include "mpif.h" - if( err.EQ.NF90_NOERR )return errmsg = NF90_STRERROR(err) print*,'' @@ -31,6 +30,8 @@ end subroutine netcdf_err subroutine error_handler(string, rc) + use mpi + implicit none character(len=*), intent(in) :: string @@ -39,8 +40,6 @@ subroutine error_handler(string, rc) integer :: ierr - include "mpif.h" - print*,"- FATAL ERROR: ", string if (present(rc)) print*,"- IOSTAT IS: ", rc call mpi_abort(mpi_comm_world, 999, ierr) From 61c25ed93ae7ac6d3e02ad3c3cfb815f82e77b52 Mon Sep 17 00:00:00 2001 From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com> Date: Thu, 1 Oct 2020 12:31:43 -0400 Subject: [PATCH 7/8] Update GDAS initialization utility for ATMOS subdirectory Update GDAS initialization utility for new 'atmos' subfolder. Add load of NDATE utility to ./gdas_init/driver.hera.sh. It was inadvertently removed during a previous merge. Issue #159. --- util/gdas_init/driver.hera.sh | 4 ++++ util/gdas_init/run_pre-v14.chgres.sh | 4 ++-- util/gdas_init/run_v14.chgres.sh | 4 ++-- util/gdas_init/run_v15.chgres.sh | 4 ++-- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/util/gdas_init/driver.hera.sh b/util/gdas_init/driver.hera.sh index 02a703e7d..3fee1147e 100755 --- a/util/gdas_init/driver.hera.sh +++ b/util/gdas_init/driver.hera.sh @@ -11,6 +11,10 @@ set -x source ../../sorc/machine-setup.sh > /dev/null 2>&1 source ../../modulefiles/build.$target +# Needed for NDATE utility +module use -a /scratch2/NCEPDEV/nwprod/NCEPLIBS/modulefiles +module load prod_util/1.1.0 + PROJECT_CODE=fv3-cpu QUEUE=batch diff --git a/util/gdas_init/run_pre-v14.chgres.sh b/util/gdas_init/run_pre-v14.chgres.sh index 92cfe80b3..42946b220 100755 --- a/util/gdas_init/run_pre-v14.chgres.sh +++ b/util/gdas_init/run_pre-v14.chgres.sh @@ -27,14 +27,14 @@ if [ ${MEMBER} == 'hires' ]; then CTAR=${CRES_HIRES} INPUT_DATA_DIR="${EXTRACT_DIR}/gdas.${yy}${mm}${dd}/${hh}" RADSTAT_DATA_DIR="${EXTRACT_DIR}/gdas.${yy}${mm}${dd}/${hh}" - OUTDIR=$OUTDIR/gdas.${yy}${mm}${dd}/${hh} + OUTDIR=$OUTDIR/gdas.${yy}${mm}${dd}/${hh}/atmos ATMFILE="gdas1.t${hh}z.sanl" SFCFILE="gdas1.t${hh}z.sfcanl" else CTAR=${CRES_ENKF} INPUT_DATA_DIR="${EXTRACT_DIR}/enkf.${yy}${mm}${dd}/${hh}/mem${MEMBER}" RADSTAT_DATA_DIR="${EXTRACT_DIR}/enkf.${yy}${mm}${dd}/${hh}/mem${MEMBER}" - OUTDIR=$OUTDIR/enkfgdas.${yy}${mm}${dd}/${hh}/mem${MEMBER} + OUTDIR=$OUTDIR/enkfgdas.${yy}${mm}${dd}/${hh}/atmos/mem${MEMBER} ATMFILE="siganl_${yy}${mm}${dd}${hh}_mem${MEMBER}" SFCFILE="sfcanl_${yy}${mm}${dd}${hh}_mem${MEMBER}" fi diff --git a/util/gdas_init/run_v14.chgres.sh b/util/gdas_init/run_v14.chgres.sh index c60c4b3d5..af35acf61 100755 --- a/util/gdas_init/run_v14.chgres.sh +++ b/util/gdas_init/run_v14.chgres.sh @@ -26,7 +26,7 @@ if [ ${MEMBER} == 'hires' ]; then CTAR=${CRES_HIRES} INPUT_DATA_DIR="${EXTRACT_DIR}/gdas.${yy}${mm}${dd}/${hh}" RADSTAT_DATA_DIR="${EXTRACT_DIR}/gdas.${yy}${mm}${dd}/${hh}" - OUTDIR=$OUTDIR/gdas.${yy}${mm}${dd}/${hh} + OUTDIR=$OUTDIR/gdas.${yy}${mm}${dd}/${hh}/atmos ATMFILE="gdas.t${hh}z.atmanl.nemsio" SFCFILE="gdas.t${hh}z.sfcanl.nemsio" NSTFILE="gdas.t${hh}z.nstanl.nemsio" @@ -34,7 +34,7 @@ else CTAR=${CRES_ENKF} INPUT_DATA_DIR="${EXTRACT_DIR}/enkf.${yy}${mm}${dd}/${hh}/mem${MEMBER}" RADSTAT_DATA_DIR="${EXTRACT_DIR}/enkf.${yy}${mm}${dd}/${hh}/mem${MEMBER}" - OUTDIR=$OUTDIR/enkfgdas.${yy}${mm}${dd}/${hh}/mem${MEMBER} + OUTDIR=$OUTDIR/enkfgdas.${yy}${mm}${dd}/${hh}/atmos/mem${MEMBER} ATMFILE="gdas.t${hh}z.ratmanl.mem${MEMBER}.nemsio" SFCFILE="gdas.t${hh}z.sfcanl.mem${MEMBER}.nemsio" NSTFILE="gdas.t${hh}z.nstanl.mem${MEMBER}.nemsio" diff --git a/util/gdas_init/run_v15.chgres.sh b/util/gdas_init/run_v15.chgres.sh index c26d92e6c..fce1569c9 100755 --- a/util/gdas_init/run_v15.chgres.sh +++ b/util/gdas_init/run_v15.chgres.sh @@ -27,13 +27,13 @@ if [ ${MEMBER} == 'hires' ]; then CTAR=${CRES_HIRES} INPUT_DATA_DIR="${EXTRACT_DIR}/gdas.${yy_d}${mm_d}${dd_d}/${hh_d}/RESTART" RADSTAT_DATA_DIR="${EXTRACT_DIR}/gdas.${yy}${mm}${dd}/${hh}" - OUTDIR=$OUTDIR/gdas.${yy}${mm}${dd}/${hh} + OUTDIR=$OUTDIR/gdas.${yy}${mm}${dd}/${hh}/atmos else CINP=C384 CTAR=${CRES_ENKF} INPUT_DATA_DIR="${EXTRACT_DIR}/enkfgdas.${yy_d}${mm_d}${dd_d}/${hh_d}/mem${MEMBER}/RESTART" RADSTAT_DATA_DIR="${EXTRACT_DIR}/enkfgdas.${yy}${mm}${dd}/${hh}/mem${MEMBER}" - OUTDIR=$OUTDIR/enkfgdas.${yy}${mm}${dd}/${hh}/mem${MEMBER} + OUTDIR=$OUTDIR/enkfgdas.${yy}${mm}${dd}/${hh}/atmos/mem${MEMBER} fi rm -fr $WORKDIR From b69f9d13f6469160f90be9d350fb60d29b3464d9 Mon Sep 17 00:00:00 2001 From: David Wright Date: Wed, 7 Oct 2020 13:48:01 -0400 Subject: [PATCH 8/8] Bug Fix for FVCOM Processing to Allow Adding NetCDF Variables on Hera Small bug fix to ./sorc/fvcom_tools.fd that allows for the addition of new variables on Hera. --- sorc/fvcom_tools.fd/module_ncio.f90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/sorc/fvcom_tools.fd/module_ncio.f90 b/sorc/fvcom_tools.fd/module_ncio.f90 index 5f118ad93..9b80d791f 100644 --- a/sorc/fvcom_tools.fd/module_ncio.f90 +++ b/sorc/fvcom_tools.fd/module_ncio.f90 @@ -2339,6 +2339,9 @@ subroutine add_new_var_3d(this,varname,dname1,dname2,dname3,lname,units) ,lname,units integer :: status, ncid, dim1id, dim2id, dim3id, varid + status = nf90_redef(this%ncid) !Enter Define Mode + if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_inq_dimid(this%ncid, dname1, dim1id) if (status /= nf90_noerr) call this%handle_err(status) status = nf90_inq_dimid(this%ncid, dname2, dim2id) @@ -2355,6 +2358,10 @@ subroutine add_new_var_3d(this,varname,dname1,dname2,dname3,lname,units) status = nf90_put_att(this%ncid, varid, 'units', units) if (status /= nf90_noerr) call this%handle_err(status) + status = nf90_enddef(this%ncid) !Exit Define Mode and + ! return to Data Mode + if (status /= nf90_noerr) call this%handle_err(status) + end subroutine add_new_var_3d end module module_ncio