diff --git a/jobs/JREGIONAL_MAKE_ICS b/jobs/JREGIONAL_MAKE_ICS index 2ddf94bac..4ca00af53 100755 --- a/jobs/JREGIONAL_MAKE_ICS +++ b/jobs/JREGIONAL_MAKE_ICS @@ -91,16 +91,6 @@ esac # #----------------------------------------------------------------------- # -# Find the directory in which the wgrib2 executable is located. -# -#----------------------------------------------------------------------- -# -wgrib2_dir=$( which wgrib2 ) || print_err_msg_exit "\ -Directory in which the wgrib2 executable is located could not be found: - wgrib2_dir = \"${wgrib2_dir}\"" -# -#----------------------------------------------------------------------- -# # Set the name of and create the directory in which the output from this # script will be placed (if that directory doesn't already exist). # @@ -117,7 +107,6 @@ mkdir_vrfy -p "${ics_dir}" # $SCRIPTSDIR/exregional_make_ics.sh \ ics_dir="${ics_dir}" \ - wgrib2_dir="${wgrib2_dir}" \ APRUN="${APRUN}" || \ print_err_msg_exit "\ Call to ex-script corresponding to J-job \"${scrfunc_fn}\" failed." diff --git a/jobs/JREGIONAL_MAKE_LBCS b/jobs/JREGIONAL_MAKE_LBCS index 49923a0db..43bd27b5d 100755 --- a/jobs/JREGIONAL_MAKE_LBCS +++ b/jobs/JREGIONAL_MAKE_LBCS @@ -92,16 +92,6 @@ esac # #----------------------------------------------------------------------- # -# Find the directory in which the wgrib2 executable is located. -# -#----------------------------------------------------------------------- -# -wgrib2_dir=$( which wgrib2 ) || print_err_msg_exit "\ -Directory in which the wgrib2 executable is located could not be found: - wgrib2_dir = \"${wgrib2_dir}\"" -# -#----------------------------------------------------------------------- -# # Set the name of and create the directory in which the output from this # script will be placed (if it doesn't already exist). # @@ -118,7 +108,6 @@ mkdir_vrfy -p "${lbcs_dir}" # $SCRIPTSDIR/exregional_make_lbcs.sh \ lbcs_dir="${lbcs_dir}" \ - wgrib2_dir="${wgrib2_dir}" \ APRUN="${APRUN}" || \ print_err_msg_exit "\ Call to ex-script corresponding to J-job \"${scrfunc_fn}\" failed." diff --git a/modulefiles/codes/cheyenne/global_equiv_resol b/modulefiles/codes/cheyenne/global_equiv_resol deleted file mode 100644 index d50a83ea0..000000000 --- a/modulefiles/codes/cheyenne/global_equiv_resol +++ /dev/null @@ -1,10 +0,0 @@ -#%Module##################################################### -## Module file for global_equiv_resol on NCAR/UCAR Cheyenne -############################################################# -module purge -module load ncarenv/1.3 -module load intel/18.0.5 -module load ncarcompilers/0.5.0 -module load netcdf/4.6.3 -# No hdf5 loaded since netcdf and hdf5 reside together on cheyenne - diff --git a/modulefiles/codes/cheyenne/mosaic_file b/modulefiles/codes/cheyenne/mosaic_file deleted file mode 100644 index a24cf1ce6..000000000 --- a/modulefiles/codes/cheyenne/mosaic_file +++ /dev/null @@ -1,10 +0,0 @@ -#%Module##################################################### -## Module file for mosaic_file on NCAR/UCAR Cheyenne -############################################################# -module purge -module load ncarenv/1.3 -module load intel/18.0.5 -module load ncarcompilers/0.5.0 -module load impi/2018.4.274 -module load netcdf/4.6.3 -# No hdf5 loaded since netcdf and hdf5 reside together on cheyenne diff --git a/modulefiles/codes/cheyenne/regional_grid b/modulefiles/codes/cheyenne/regional_grid deleted file mode 100644 index bd3d7874c..000000000 --- a/modulefiles/codes/cheyenne/regional_grid +++ /dev/null @@ -1,10 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid on NCAR/UCAR Cheyenne -############################################################# -module purge -module load ncarenv/1.3 -module load intel/18.0.5 -module load ncarcompilers/0.5.0 -module load impi/2018.4.274 -module load netcdf/4.6.3 -# No hdf5 loaded since netcdf and hdf5 reside together on cheyenne diff --git a/modulefiles/codes/hera/global_equiv_resol b/modulefiles/codes/hera/global_equiv_resol deleted file mode 100644 index 532e7fb65..000000000 --- a/modulefiles/codes/hera/global_equiv_resol +++ /dev/null @@ -1,7 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module purge -module load intel/18.0.5.274 -module load netcdf/4.6.1 -module load hdf5/1.10.4 diff --git a/modulefiles/codes/hera/mosaic_file b/modulefiles/codes/hera/mosaic_file deleted file mode 100644 index 532e7fb65..000000000 --- a/modulefiles/codes/hera/mosaic_file +++ /dev/null @@ -1,7 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module purge -module load intel/18.0.5.274 -module load netcdf/4.6.1 -module load hdf5/1.10.4 diff --git a/modulefiles/codes/hera/regional_grid b/modulefiles/codes/hera/regional_grid deleted file mode 100644 index 532e7fb65..000000000 --- a/modulefiles/codes/hera/regional_grid +++ /dev/null @@ -1,7 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module purge -module load intel/18.0.5.274 -module load netcdf/4.6.1 -module load hdf5/1.10.4 diff --git a/modulefiles/codes/jet/global_equiv_resol b/modulefiles/codes/jet/global_equiv_resol deleted file mode 100644 index 532e7fb65..000000000 --- a/modulefiles/codes/jet/global_equiv_resol +++ /dev/null @@ -1,7 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module purge -module load intel/18.0.5.274 -module load netcdf/4.6.1 -module load hdf5/1.10.4 diff --git a/modulefiles/codes/jet/mosaic_file b/modulefiles/codes/jet/mosaic_file deleted file mode 100644 index 532e7fb65..000000000 --- a/modulefiles/codes/jet/mosaic_file +++ /dev/null @@ -1,7 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module purge -module load intel/18.0.5.274 -module load netcdf/4.6.1 -module load hdf5/1.10.4 diff --git a/modulefiles/codes/jet/regional_grid b/modulefiles/codes/jet/regional_grid deleted file mode 100644 index c5efa3218..000000000 --- a/modulefiles/codes/jet/regional_grid +++ /dev/null @@ -1,8 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module purge -module load intel/18.0.5.274 -module load netcdf/4.7.0 -module load hdf5/1.10.5 - diff --git a/modulefiles/codes/odin/global_equiv_resol b/modulefiles/codes/odin/global_equiv_resol deleted file mode 100644 index 123d202ea..000000000 --- a/modulefiles/codes/odin/global_equiv_resol +++ /dev/null @@ -1,13 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module load PrgEnv-intel -module swap intel/19.0.5.281 -module load cray-mpich/7.7.10 -module load cray-libsci -module load cray-netcdf-hdf5parallel -module load cray-parallel-netcdf -module load cray-hdf5-parallel - -export NETCDF=/opt/cray/pe/netcdf-hdf5parallel/4.6.3.2/INTEL/19.0 -export HDF5=/opt/cray/pe/hdf5-parallel/1.10.5.2/INTEL/19.0 diff --git a/modulefiles/codes/odin/mosaic_file b/modulefiles/codes/odin/mosaic_file deleted file mode 100644 index 123d202ea..000000000 --- a/modulefiles/codes/odin/mosaic_file +++ /dev/null @@ -1,13 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module load PrgEnv-intel -module swap intel/19.0.5.281 -module load cray-mpich/7.7.10 -module load cray-libsci -module load cray-netcdf-hdf5parallel -module load cray-parallel-netcdf -module load cray-hdf5-parallel - -export NETCDF=/opt/cray/pe/netcdf-hdf5parallel/4.6.3.2/INTEL/19.0 -export HDF5=/opt/cray/pe/hdf5-parallel/1.10.5.2/INTEL/19.0 diff --git a/modulefiles/codes/odin/regional_grid b/modulefiles/codes/odin/regional_grid deleted file mode 100644 index 2d4b9a9a3..000000000 --- a/modulefiles/codes/odin/regional_grid +++ /dev/null @@ -1,11 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module load PrgEnv-intel -module swap intel/19.0.5.281 -module load cray-mpich/7.7.10 -module load cray-libsci -module load cray-netcdf-hdf5parallel -module load cray-parallel-netcdf -module load cray-hdf5-parallel - diff --git a/modulefiles/codes/stampede/global_equiv_resol b/modulefiles/codes/stampede/global_equiv_resol deleted file mode 100644 index 0f78231fa..000000000 --- a/modulefiles/codes/stampede/global_equiv_resol +++ /dev/null @@ -1,10 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module load intel/18.0.2 -module load impi/18.0.2 -module load hdf5/1.10.4 -module load netcdf/4.6.2 - -setenv NETCDF /opt/apps/intel18/netcdf/4.6.2/x86_64 -setenv HDF5 /opt/apps/intel18/hdf5/1.10.4/x86_64 diff --git a/modulefiles/codes/stampede/mosaic_file b/modulefiles/codes/stampede/mosaic_file deleted file mode 100644 index bab633db7..000000000 --- a/modulefiles/codes/stampede/mosaic_file +++ /dev/null @@ -1,12 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# - - -module load intel/18.0.2 -module load impi/18.0.2 -module load hdf5/1.10.4 -module load netcdf/4.6.2 - -setenv NETCDF /opt/apps/intel18/netcdf/4.6.2/x86_64 -setenv HDF5 /opt/apps/intel18/hdf5/1.10.4/x86_64 diff --git a/modulefiles/codes/stampede/regional_grid b/modulefiles/codes/stampede/regional_grid deleted file mode 100644 index 0f78231fa..000000000 --- a/modulefiles/codes/stampede/regional_grid +++ /dev/null @@ -1,10 +0,0 @@ -#%Module##################################################### -## Module file for regional_grid -############################################################# -module load intel/18.0.2 -module load impi/18.0.2 -module load hdf5/1.10.4 -module load netcdf/4.6.2 - -setenv NETCDF /opt/apps/intel18/netcdf/4.6.2/x86_64 -setenv HDF5 /opt/apps/intel18/hdf5/1.10.4/x86_64 diff --git a/modulefiles/tasks/cheyenne/make_grid b/modulefiles/tasks/cheyenne/make_grid deleted file mode 100644 index aeb62bab7..000000000 --- a/modulefiles/tasks/cheyenne/make_grid +++ /dev/null @@ -1,15 +0,0 @@ -#%Module##################################################### - -module purge - -module load ncarenv/1.3 -module load intel/19.0.2 -module load ncarcompilers/0.5.0 -module load netcdf/4.6.3 - -if [module-info mode load] { - system "ncar_pylib /glade/p/ral/jntp/UFS_CAM/ncar_pylib_20200427" -} -if [module-info mode remove] { - system "deactivate" -} diff --git a/modulefiles/tasks/cheyenne/make_grid.local b/modulefiles/tasks/cheyenne/make_grid.local new file mode 100644 index 000000000..7a83a9d8d --- /dev/null +++ b/modulefiles/tasks/cheyenne/make_grid.local @@ -0,0 +1,9 @@ + +if [module-info mode load] { + system "ncar_pylib /glade/p/ral/jntp/UFS_CAM/ncar_pylib_20200427" +} + +if [module-info mode remove] { + system "deactivate" +} + diff --git a/modulefiles/tasks/cheyenne/make_ics.local b/modulefiles/tasks/cheyenne/make_ics.local index 4a271536c..7a83a9d8d 100644 --- a/modulefiles/tasks/cheyenne/make_ics.local +++ b/modulefiles/tasks/cheyenne/make_ics.local @@ -1,2 +1,9 @@ -ncar_pylib /glade/p/ral/jntp/UFS_CAM/ncar_pylib_20200427 +if [module-info mode load] { + system "ncar_pylib /glade/p/ral/jntp/UFS_CAM/ncar_pylib_20200427" +} + +if [module-info mode remove] { + system "deactivate" +} + diff --git a/modulefiles/tasks/cheyenne/make_lbcs.local b/modulefiles/tasks/cheyenne/make_lbcs.local index 4a271536c..7a83a9d8d 100644 --- a/modulefiles/tasks/cheyenne/make_lbcs.local +++ b/modulefiles/tasks/cheyenne/make_lbcs.local @@ -1,2 +1,9 @@ -ncar_pylib /glade/p/ral/jntp/UFS_CAM/ncar_pylib_20200427 +if [module-info mode load] { + system "ncar_pylib /glade/p/ral/jntp/UFS_CAM/ncar_pylib_20200427" +} + +if [module-info mode remove] { + system "deactivate" +} + diff --git a/modulefiles/tasks/cheyenne/run_fcst.local b/modulefiles/tasks/cheyenne/run_fcst.local index 5dcc46add..7a83a9d8d 100644 --- a/modulefiles/tasks/cheyenne/run_fcst.local +++ b/modulefiles/tasks/cheyenne/run_fcst.local @@ -1,2 +1,9 @@ -system "ncar_pylib /glade/p/ral/jntp/UFS_CAM/ncar_pylib_20200427" +if [module-info mode load] { + system "ncar_pylib /glade/p/ral/jntp/UFS_CAM/ncar_pylib_20200427" +} + +if [module-info mode remove] { + system "deactivate" +} + diff --git a/modulefiles/tasks/hera/make_grid b/modulefiles/tasks/hera/make_grid deleted file mode 100644 index 12dc7ae10..000000000 --- a/modulefiles/tasks/hera/make_grid +++ /dev/null @@ -1,15 +0,0 @@ -#%Module##################################################### -## Module file for make_grid task. -############################################################# - -module purge - -module load intel/18.0.5.274 -module load netcdf/4.7.0 -module load hdf5/1.10.5 - -module use -a /contrib/miniconda3/modulefiles -module load miniconda3 -if [module-info mode load] { - system "conda activate regional_workflow" -} diff --git a/modulefiles/tasks/hera/make_grid.local b/modulefiles/tasks/hera/make_grid.local new file mode 100644 index 000000000..e7d11ba40 --- /dev/null +++ b/modulefiles/tasks/hera/make_grid.local @@ -0,0 +1,6 @@ + +module use -a /contrib/miniconda3/modulefiles +module load miniconda3 +if [module-info mode load] { + system "conda activate regional_workflow" +} diff --git a/modulefiles/tasks/hera/make_ics.local b/modulefiles/tasks/hera/make_ics.local index ff1cad391..e7d11ba40 100644 --- a/modulefiles/tasks/hera/make_ics.local +++ b/modulefiles/tasks/hera/make_ics.local @@ -1,5 +1,6 @@ -module load wgrib2 module use -a /contrib/miniconda3/modulefiles module load miniconda3 -conda activate regional_workflow +if [module-info mode load] { + system "conda activate regional_workflow" +} diff --git a/modulefiles/tasks/hera/make_lbcs.local b/modulefiles/tasks/hera/make_lbcs.local index ff1cad391..e7d11ba40 100644 --- a/modulefiles/tasks/hera/make_lbcs.local +++ b/modulefiles/tasks/hera/make_lbcs.local @@ -1,5 +1,6 @@ -module load wgrib2 module use -a /contrib/miniconda3/modulefiles module load miniconda3 -conda activate regional_workflow +if [module-info mode load] { + system "conda activate regional_workflow" +} diff --git a/modulefiles/tasks/jet/make_grid b/modulefiles/tasks/jet/make_grid deleted file mode 100644 index 12dc7ae10..000000000 --- a/modulefiles/tasks/jet/make_grid +++ /dev/null @@ -1,15 +0,0 @@ -#%Module##################################################### -## Module file for make_grid task. -############################################################# - -module purge - -module load intel/18.0.5.274 -module load netcdf/4.7.0 -module load hdf5/1.10.5 - -module use -a /contrib/miniconda3/modulefiles -module load miniconda3 -if [module-info mode load] { - system "conda activate regional_workflow" -} diff --git a/modulefiles/tasks/jet/make_grid.local b/modulefiles/tasks/jet/make_grid.local new file mode 100644 index 000000000..e7d11ba40 --- /dev/null +++ b/modulefiles/tasks/jet/make_grid.local @@ -0,0 +1,6 @@ + +module use -a /contrib/miniconda3/modulefiles +module load miniconda3 +if [module-info mode load] { + system "conda activate regional_workflow" +} diff --git a/modulefiles/tasks/jet/make_ics.local b/modulefiles/tasks/jet/make_ics.local index d5fa8e56b..d5e37266b 100644 --- a/modulefiles/tasks/jet/make_ics.local +++ b/modulefiles/tasks/jet/make_ics.local @@ -2,4 +2,6 @@ module load wgrib2/2.0.8 module use -a /contrib/miniconda3/modulefiles module load miniconda3 -conda activate regional_workflow +if [module-info mode load] { + system "conda activate regional_workflow" +} diff --git a/modulefiles/tasks/jet/make_lbcs.local b/modulefiles/tasks/jet/make_lbcs.local index d5fa8e56b..d5e37266b 100644 --- a/modulefiles/tasks/jet/make_lbcs.local +++ b/modulefiles/tasks/jet/make_lbcs.local @@ -2,4 +2,6 @@ module load wgrib2/2.0.8 module use -a /contrib/miniconda3/modulefiles module load miniconda3 -conda activate regional_workflow +if [module-info mode load] { + system "conda activate regional_workflow" +} diff --git a/modulefiles/tasks/odin/make_grid b/modulefiles/tasks/odin/make_grid deleted file mode 100644 index d37cb0e9d..000000000 --- a/modulefiles/tasks/odin/make_grid +++ /dev/null @@ -1,14 +0,0 @@ -#%Module##################################################### -## Module file for make_grid task. -############################################################# - -#module purge - -module load PrgEnv-intel -module swap intel/19.0.5.281 -module load cray-libsci -module load cray-netcdf-hdf5parallel -module load cray-parallel-netcdf -module load cray-hdf5-parallel -module load gcc/6.1.0 -module load slurm diff --git a/modulefiles/tasks/stampede/make_grid b/modulefiles/tasks/stampede/make_grid deleted file mode 100644 index aa867cfbf..000000000 --- a/modulefiles/tasks/stampede/make_grid +++ /dev/null @@ -1,11 +0,0 @@ -#%Module##################################################### -## Module file for make_grid task. -############################################################# - -module purge - -module load intel/18.0.2 -module load impi/18.0.2 -module load hdf5/1.10.4 -module load netcdf/4.6.2 -module load python3/3.7.0 diff --git a/modulefiles/tasks/wcoss_cray/make_grid.local b/modulefiles/tasks/wcoss_cray/make_grid.local new file mode 100644 index 000000000..4bebfc5e7 --- /dev/null +++ b/modulefiles/tasks/wcoss_cray/make_grid.local @@ -0,0 +1,33 @@ +module load modules + +module load xt-lsfhpc +module load ncep +module load alps +module load dvs +module load xpmem +module load ugni +module load craype-network-aries +module load switch +module load rca +module load gni-headers +module load udreg +module load hpss + +module load prod_util +module load g2tmpl-intel/1.4.0 +module load crtm-intel/2.2.6 +module load iobuf/2.0.7 +module load gempak/7.3.0 + +module load nco-gnu-sandybridge/4.4.4 +module load NetCDF-intel-sandybridge/4.2 +module load cfp-intel-sandybridge/1.1.0 +export USE_CFP=YES + +module load grib_util/1.1.0 + +module use -a /gpfs/hps3/emc/nems/noscrub/emc.nemspara/soft/modulefiles +module load esmf/8.0.0 + +module load python/3.6.3 + diff --git a/modulefiles/tasks/wcoss_dell_p3/make_grid.local b/modulefiles/tasks/wcoss_dell_p3/make_grid.local new file mode 100644 index 000000000..b93479d8c --- /dev/null +++ b/modulefiles/tasks/wcoss_dell_p3/make_grid.local @@ -0,0 +1,2 @@ +module load lsf/10.1 +module load python/3.6.3 diff --git a/scripts/exregional_make_grid.sh b/scripts/exregional_make_grid.sh index 977c9c99c..2128e9ddb 100755 --- a/scripts/exregional_make_grid.sh +++ b/scripts/exregional_make_grid.sh @@ -372,9 +372,7 @@ generation executable (exec_fp): 'dely': ${DEL_ANGLE_Y_SG}, 'lx': ${NEG_NX_OF_DOM_WITH_WIDE_HALO}, 'ly': ${NEG_NY_OF_DOM_WITH_WIDE_HALO}, - 'a': ${ESGgrid_ALPHA_PARAM}, - 'k': ${ESGgrid_KAPPA_PARAM}, -} + } " # # Call the python script to create the namelist file. diff --git a/scripts/exregional_make_ics.sh b/scripts/exregional_make_ics.sh index 661918690..6759edc64 100755 --- a/scripts/exregional_make_ics.sh +++ b/scripts/exregional_make_ics.sh @@ -44,7 +44,7 @@ In directory: \"${scrfunc_dir}\" This is the ex-script for the task that generates initial condition (IC), surface, and zeroth hour lateral boundary condition (LBC0) files -for FV3 (in NetCDF format). +(in NetCDF format) for the FV3-LAM. ========================================================================" # #----------------------------------------------------------------------- @@ -56,7 +56,6 @@ for FV3 (in NetCDF format). #----------------------------------------------------------------------- # valid_args=( \ -"wgrib2_dir" \ "ics_dir" \ "APRUN" \ ) @@ -95,34 +94,31 @@ cd_vrfy $workdir # #----------------------------------------------------------------------- # -# Set physics-suite-dependent variables that are needed in the FORTRAN -# namelist file that the chgres executable will read in. +# Set physics-suite-dependent variable mapping table needed in the FORTRAN +# namelist file that the chgres_cube executable will read in. # #----------------------------------------------------------------------- # -phys_suite="" +varmap_file="" case "${CCPP_PHYS_SUITE}" in -"FV3_GFS_2017_gfdlmp" | "FV3_GFS_2017_gfdlmp_regional" ) - phys_suite="GFS" +"FV3_GFS_2017_gfdlmp" | "FV3_GFS_2017_gfdlmp_regional" | "FV3_GFS_v16beta" | \ +"FV3_GFS_v15p2" | "FV3_CPT_v0" ) + varmap_file="GFSphys_var_map.txt" ;; -"FV3_GSD_v0" | "FV3_GSD_SAR" | "FV3_RRFS_v1beta") - phys_suite="GSD" - ;; -"FV3_CPT_v0") - phys_suite="CPT" - ;; -"FV3_GFS_v15p2") - phys_suite="v15p2" - ;; -"FV3_GFS_v16beta") - phys_suite="v16beta" +"FV3_GSD_v0" | "FV3_GSD_SAR" | \ +"FV3_RRFS_v1beta" ) + if [ "${EXTRN_MDL_NAME_ICS}" = "RAPX" ] || [ "${EXTRN_MDL_NAME_ICS}" = "HRRRX" ]; then + varmap_file="GSDphys_var_map.txt" + elif [ "${EXTRN_MDL_NAME_ICS}" = "NAM" ] || [ "${EXTRN_MDL_NAME_ICS}" = "FV3GFS" ] || \ + [ "${EXTRN_MDL_NAME_ICS}" = "GSMGFS" ]; then + varmap_file="GFSphys_var_map.txt" + fi ;; *) print_err_msg_exit "\ -Physics-suite-dependent namelist variables have not yet been specified -for this physics suite: +A variable mapping table has not yet been defined for this physics suite: CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" ;; @@ -131,7 +127,7 @@ esac #----------------------------------------------------------------------- # # Set external-model-dependent variables that are needed in the FORTRAN -# namelist file that the chgres executable will read in. These are de- +# namelist file that the chgres_cube executable will read in. These are de- # scribed below. Note that for a given external model, usually only a # subset of these all variables are set (since some may be irrelevant). # @@ -148,10 +144,10 @@ esac # model that contains the surface fields. # # input_type: -# The "type" of input being provided to chgres. This contains a combi- +# The "type" of input being provided to chgres_cube. This contains a combi- # nation of information on the external model, external model file for- # mat, and maybe other parameters. For clarity, it would be best to -# eliminate this variable in chgres and replace with with 2 or 3 others +# eliminate this variable in chgres_cube and replace with with 2 or 3 others # (e.g. extrn_mdl, extrn_mdl_file_format, etc). # # tracers_input: @@ -167,26 +163,24 @@ esac # ment of tracers should be the name to use for the O3 mixing ratio in # the output file. For GSD physics, three additional tracers -- ice, # rain, and water number concentrations -- may be specified at the end -# of tracers, and these will be calculated by chgres. -# -# internal_GSD: -# Logical variable indicating whether or not to try to read in land sur- -# face model (LSM) variables available in the HRRRX grib2 files created -# after about 2019111500. +# of tracers, and these will be calculated by chgres_cube. # -# numsoil_out: +# nsoill_out: # The number of soil layers to include in the output NetCDF file. # -# replace_FIELD, where FIELD="vgtyp", "sotyp", or "vgfrc": -# Logical variable indicating whether or not to obtain the field in +# FIELD_from_climo, where FIELD = "vgtyp", "sotyp", "vgfrc", "lai", or +# "minmax_vgfrc": +# Logical variable indicating whether or not to obtain the field in # question from climatology instead of the external model. The field in # question is one of vegetation type (FIELD="vgtyp"), soil type (FIELD= -# "sotyp"), and vegetation fraction (FIELD="vgfrc"). If replace_FIELD -# is set to ".true.", then the field is obtained from climatology (re- -# gardless of whether or not it exists in an external model file). If -# it is set to ".false.", then the field is obtained from the external -# model. If the external model file does not provide this field, then -# chgres prints out an error message and stops. +# "sotyp"), vegetation fraction (FIELD="vgfrc"), leaf area index +# (FIELD="lai"), or min/max areal fractional coverage of annual green +# vegetation (FIELD="minmax_vfrr"). If FIELD_from_climo is set to +# ".true.", then the field is obtained from climatology (regardless of +# whether or not it exists in an external model file). If it is set +# to ".false.", then the field is obtained from the external model. +# If "false" is chosen and the external model file does not provide +# this field, then chgres_cube prints out an error message and stops. # # tg3_from_soil: # Logical variable indicating whether or not to set the tg3 soil tempe- # Needs to be verified. @@ -228,7 +222,6 @@ esac # A non-prognostic variable that appears in the field_table for GSD physics # is cld_amt. Why is that in the field_table at all (since it is a non- # prognostic field), and how should we handle it here?? - # I guess this works for FV3GFS but not for the spectral GFS since these # variables won't exist in the spectral GFS atmanl files. # tracers_input="\"sphum\",\"liq_wat\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\",\"o3mr\"" @@ -243,15 +236,16 @@ fn_grib2="" input_type="" tracers_input="\"\"" tracers="\"\"" -internal_GSD="" -numsoil_out="" +nsoill_out="" geogrid_file_input_grid="\"\"" -replace_vgtyp="" -replace_sotyp="" -replace_vgfrc="" +vgtyp_from_climo="" +sotyp_from_climo="" +vgfrc_from_climo="" +minmax_vgfrc_from_climo="" +lai_from_climo="" tg3_from_soil="" convert_nst="" - +thomp_mp_climo_file="" case "${EXTRN_MDL_NAME_ICS}" in @@ -262,70 +256,73 @@ case "${EXTRN_MDL_NAME_ICS}" in fn_atm_nemsio="${EXTRN_MDL_FNS[0]}" fn_sfc_nemsio="${EXTRN_MDL_FNS[1]}" - input_type="gfs_gaussian" # For spectral GFS Gaussian grid in nemsio format. + input_type="gfs_gaussian_nemsio" # For spectral GFS Gaussian grid in nemsio format. + convert_nst=False tracers_input="[\"spfh\",\"clwmr\",\"o3mr\"]" tracers="[\"sphum\",\"liq_wat\",\"o3mr\"]" +# +# Use Thompson climatology for ice- and water-friendly aerosols if CCPP suite uses Thompson MP +# + if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then + thomp_mp_climo_file="" + elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then + thomp_mp_climo_file="${FIXam}/Thompson_MP_MONTHLY_CLIMO.nc" + else + print_err_msg_exit "\ + The chosen CCPP physics suite is unsupported at this time: + CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" + fi - internal_GSD=False - numsoil_out="4" - replace_vgtyp=True - replace_sotyp=True - replace_vgfrc=True + nsoill_out="4" #If the CCPP suites uses RUC-LSM, the scheme will interpolate from 4 to 9 soil levels. + vgtyp_from_climo=True + sotyp_from_climo=True + vgfrc_from_climo=True + minmax_vgfrc_from_climo=True + lai_from_climo=True tg3_from_soil=False - convert_nst=False ;; "FV3GFS") +# +# Use Thompson climatology for ice- and water-friendly aerosols if CCPP suite uses Thompson MP +# + if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then + thomp_mp_climo_file="" + elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then + thomp_mp_climo_file="${FIXam}/Thompson_MP_MONTHLY_CLIMO.nc" + else + print_err_msg_exit "\ + The chosen CCPP physics suite is unsupported as this time: + CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" + fi + if [ "${FV3GFS_FILE_FMT_ICS}" = "nemsio" ]; then external_model="FV3GFS" - fn_atm_nemsio="${EXTRN_MDL_FNS[0]}" - fn_sfc_nemsio="${EXTRN_MDL_FNS[1]}" - input_type="gaussian" # For FV3-GFS Gaussian grid in nemsio format. - tracers_input="[\"spfh\",\"clwmr\",\"o3mr\",\"icmr\",\"rwmr\",\"snmr\",\"grle\"]" + tracers="[\"sphum\",\"liq_wat\",\"o3mr\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\"]" -# -# If CCPP is being used, then the list of atmospheric tracers to include -# in the output file depends on the physics suite. Hopefully, this me- -# thod of specifying output tracers will be replaced with a variable -# table (which should be specific to each combination of external model, -# external model file type, and physics suite). -# - if [ "${USE_CCPP}" = "TRUE" ]; then - if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then - tracers="[\"sphum\",\"liq_wat\",\"o3mr\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\"]" - elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then -# For GSD physics, add three additional tracers (the ice, rain and water -# number concentrations) that are required for Thompson microphysics. - tracers="[\"sphum\",\"liq_wat\",\"o3mr\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\",\"ice_nc\",\"rain_nc\",\"water_nc\"]" - else - print_err_msg_exit "\ -The parameter \"tracers\" has not been defined for this combination of -external model (EXTRN_MDL_NAME_ICS), physics suite (CCPP_PHYS_SUITE), and -FV3GFS file type (FV3GFS_FILE_FMT_ICS): - EXTRN_MDL_NAME_ICS = \"${EXTRN_MDL_NAME_ICS}\" - CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\" - FV3GFS_FILE_FMT_ICS = \"${FV3GFS_FILE_FMT_ICS}\"" - fi -# -# If CCPP is not being used, the only physics suite that can be used is -# GFS. -# - else - tracers="[\"sphum\",\"liq_wat\",\"o3mr\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\"]" - fi + fn_atm_nemsio="${EXTRN_MDL_FNS[0]}" + fn_sfc_nemsio="${EXTRN_MDL_FNS[1]}" + input_type="gaussian_nemsio" # For FV3-GFS Gaussian grid in nemsio format. + convert_nst=True elif [ "${FV3GFS_FILE_FMT_ICS}" = "grib2" ]; then @@ -333,16 +330,17 @@ FV3GFS file type (FV3GFS_FILE_FMT_ICS): fn_grib2="${EXTRN_MDL_FNS[0]}" input_type="grib2" - + convert_nst=False + fi - - internal_GSD=False - numsoil_out="4" - replace_vgtyp=True - replace_sotyp=True - replace_vgfrc=True + + nsoill_out="4" #If the CCPP suites uses RUC-LSM, the scheme will interpolate from 4 to 9 soil levels. + vgtyp_from_climo=True + sotyp_from_climo=True + vgfrc_from_climo=True + minmax_vgfrc_from_climo=True + lai_from_climo=True tg3_from_soil=False - convert_nst=True ;; @@ -353,36 +351,21 @@ FV3GFS file type (FV3GFS_FILE_FMT_ICS): fn_grib2="${EXTRN_MDL_FNS[0]}" input_type="grib2" - - internal_GSD=False - cdate_min_HRRRX="2019111500" - if [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" -o \ - "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ] && \ - [ ${CDATE} -gt ${cdate_min_HRRRX} ]; then - print_info_msg " -Setting the chgres_cube namelist setting \"internal_GSD\" to \".true.\" in -order to read in land surface model (LSM) variables available in the -HRRRX grib2 files created after about \"${cdate_min_HRRRX}\"..." - internal_GSD=True - fi - +# +# Set soil levels based on LSM in CCPP SDF (RUC-LSM or Noah/Noah MP) +# if [ "${USE_CCPP}" = "TRUE" ]; then if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then - numsoil_out="4" + nsoill_out="4" elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then - numsoil_out="9" - else - print_err_msg_exit "\ -The parameter \"numsoil_out\" has not been defined for this combination -of external model (EXTRN_MDL_NAME_ICS) and physics suite (CCPP_PHYS_SUITE): - EXTRN_MDL_NAME_ICS = \"${EXTRN_MDL_NAME_ICS}\" - CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" + nsoill_out="9" fi fi # @@ -392,13 +375,16 @@ of external model (EXTRN_MDL_NAME_ICS) and physics suite (CCPP_PHYS_SUITE): geogrid_file_input_grid="/scratch2/BMC/det/beck/FV3-LAM/geo_em.d01.nc_HRRRX" elif [ "${MACHINE}" = "JET" ]; then geogrid_file_input_grid="/misc/whome/rtrr/HRRR/static/WPS/geo_em.d01.nc" - elif [ "${MACHINE}" = "CHEYENNE" ]; then - geogrid_file_input_grid="/glade/p/ral/jntp/UFS_CAM/fix/geo_em.d01.nc_HRRRX" fi - replace_vgtyp=False - replace_sotyp=False - replace_vgfrc=False + #Note that vgfrc, shdmin/shdmax (minmax_vgfrc), and lai fields are only available in HRRRX + #files after mid-July 2019, and only so long as the record order didn't change afterward + + vgtyp_from_climo=True + sotyp_from_climo=True + vgfrc_from_climo=True + minmax_vgfrc_from_climo=True + lai_from_climo=True tg3_from_soil=True convert_nst=False @@ -410,9 +396,9 @@ of external model (EXTRN_MDL_NAME_ICS) and physics suite (CCPP_PHYS_SUITE): fn_grib2="${EXTRN_MDL_FNS[0]}" input_type="grib2" - - internal_GSD=False - +# +# Set soil levels based on CCPP SDF +# if [ "${USE_CCPP}" = "TRUE" ]; then if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ @@ -420,16 +406,10 @@ of external model (EXTRN_MDL_NAME_ICS) and physics suite (CCPP_PHYS_SUITE): [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then - numsoil_out="4" + nsoill_out="4" elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then - numsoil_out="9" - else - print_err_msg_exit "\ -The parameter \"numsoil_out\" has not been defined for this combination -of external model (EXTRN_MDL_NAME_ICS) and physics suite (CCPP_PHYS_SUITE): - EXTRN_MDL_NAME_ICS = \"${EXTRN_MDL_NAME_ICS}\" - CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" + nsoill_out="9" fi fi # @@ -443,14 +423,51 @@ of external model (EXTRN_MDL_NAME_ICS) and physics suite (CCPP_PHYS_SUITE): geogrid_file_input_grid="/glade/p/ral/jntp/UFS_CAM/fix/geo_em.d01.nc_RAPX" fi - replace_vgtyp=False - replace_sotyp=False - replace_vgfrc=False + vgtyp_from_climo=True + sotyp_from_climo=False + vgfrc_from_climo=True + minmax_vgfrc_from_climo=True + lai_from_climo=True tg3_from_soil=True convert_nst=False ;; +"NAM") + + external_model="NAM" + + fn_grib2="${EXTRN_MDL_FNS[0]}" + input_type="grib2" + +# +# Use Thompson climatology for ice- and water-friendly aerosols if CCPP suite uses Thompson MP +# + if [ "${USE_CCPP}" = "TRUE" ]; then + if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then + thomp_mp_climo_file="" + elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then + [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ + thomp_mp_climo_file="${FIXam}/Thompson_MP_MONTHLY_CLIMO.nc" + fi + fi + + nsoill_out="4" #If the CCPP suites uses RUC-LSM, the scheme will interpolate from 4 to 9 soil levels. + vgtyp_from_climo=True + sotyp_from_climo=True + vgfrc_from_climo=True + minmax_vgfrc_from_climo=True + lai_from_climo=True + tg3_from_soil=False + convert_nst=False + + ;; + *) print_err_msg_exit "\ External-model-dependent namelist variables have not yet been specified @@ -491,49 +508,6 @@ fi # Build the FORTRAN namelist file that chgres_cube will read in. # #----------------------------------------------------------------------- -# -# For GFS physics, the character arrays tracers_input(:) and tracers(:) -# must be specified in the namelist file. tracers_input(:) contains the -# tracer name to look for in the external model file(s), while tracers(:) -# contains the names to use for the tracers in the output NetCDF files -# that chgres creates (that will be read in by FV3). Since when FV3 -# reads these NetCDF files it looks for atmospheric traces as specified -# in the file field_table, tracers(:) should be set to the names in -# field_table. -# -# NOTE: This process should be automated where the set of elements that -# tracers(:) should be set to is obtained from reading in field_table. -# -# To know how to set tracers_input(:), you have to know the names of the -# variables in the input atmospheric nemsio file (usually this file is -# named gfs.t00z.atmanl.nemsio). -# -# It is not quite clear how these should be specified. Here are a list -# of examples: -# -# [Gerard.Ketefian@tfe05] /scratch3/.../chgres_cube.fd/run (feature/chgres_grib2_gsk) -# $ grep -n -i "tracers" * | grep theia -# config.C1152.l91.atm.theia.nml:24: tracers="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" -# config.C1152.l91.atm.theia.nml:25: tracers_input="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" -# config.C48.gaussian.theia.nml:20: tracers="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" -# config.C48.gaussian.theia.nml:21: tracers_input="spfh","clwmr","o3mr","icmr","rwmr","snmr","grle" -# config.C48.gfs.gaussian.theia.nml:21: tracers="sphum","liq_wat","o3mr" -# config.C48.gfs.gaussian.theia.nml:22: tracers_input="spfh","clwmr","o3mr" -# config.C48.gfs.spectral.theia.nml:21: tracers_input="spfh","o3mr","clwmr" -# config.C48.gfs.spectral.theia.nml:22: tracers="sphum","o3mr","liq_wat" -# config.C48.theia.nml:21: tracers="sphum","liq_wat","o3mr" -# config.C48.theia.nml:22: tracers_input="spfh","clwmr","o3mr" -# config.C768.atm.theia.nml:24: tracers="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" -# config.C768.atm.theia.nml:25: tracers_input="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" -# config.C768.l91.atm.theia.nml:24: tracers="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" -# config.C768.l91.atm.theia.nml:25: tracers_input="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" -# config.C768.nest.atm.theia.nml:22: tracers="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" -# config.C768.nest.atm.theia.nml:23: tracers_input="sphum","liq_wat","o3mr","ice_wat","rainwat","snowwat","graupel" - - -# fix_dir_target_grid="${BASEDIR}/ESG_grid_HRRR_like_fix_files_chgres_cube" -# base_install_dir="${SORCDIR}/chgres_cube.fd" - # # Create a multiline variable that consists of a yaml-compliant string # specifying the values that the namelist variables need to be set to @@ -543,39 +517,39 @@ fi # settings=" 'config': { + 'fix_dir_input_grid': ${FIXgsm}, 'fix_dir_target_grid': ${FIXLAM}, - 'mosaic_file_target_grid': ${FIXLAM}/${CRES}${DOT_OR_USCORE}mosaic.halo${NH4}.nc, + 'mosaic_file_target_grid': ${FIXLAM}/${CRES}${DOT_OR_USCORE}mosaic.halo$((10#${NH4})).nc, 'orog_dir_target_grid': ${FIXLAM}, - 'orog_files_target_grid': ${CRES}${DOT_OR_USCORE}oro_data.tile${TILE_RGNL}.halo${NH4}.nc, + 'orog_files_target_grid': ${CRES}${DOT_OR_USCORE}oro_data.tile${TILE_RGNL}.halo$((10#${NH4})).nc, 'vcoord_file_target_grid': ${FIXam}/global_hyblev.l65.txt, - 'mosaic_file_input_grid': '', - 'orog_dir_input_grid': '', - 'base_install_dir': ${CHGRES_DIR}, - 'wgrib2_path': ${wgrib2_dir}, + 'varmap_file': ${UFS_UTILS_DIR}/parm/varmap_tables/${varmap_file}, 'data_dir_input_grid': ${extrn_mdl_staging_dir}, 'atm_files_input_grid': ${fn_atm_nemsio}, 'sfc_files_input_grid': ${fn_sfc_nemsio}, 'grib2_file_input_grid': \"${fn_grib2}\", - 'cycle_mon': $((10#$mm)), - 'cycle_day': $((10#$dd)), - 'cycle_hour': $((10#$hh)), + 'cycle_mon': $((10#${mm})), + 'cycle_day': $((10#${dd})), + 'cycle_hour': $((10#${hh})), 'convert_atm': True, 'convert_sfc': True, 'convert_nst': ${convert_nst}, 'regional': 1, - 'halo_bndy': ${NH4}, + 'halo_bndy': $((10#${NH4})), + 'halo_blend': $((10#${HALO_BLEND})), 'input_type': ${input_type}, 'external_model': ${external_model}, 'tracers_input': ${tracers_input}, - 'tracers': ${tracers}, - 'phys_suite': ${phys_suite}, - 'internal_GSD': ${internal_GSD}, - 'numsoil_out': ${numsoil_out}, + 'tracers': ${tracers}, + 'nsoill_out': $((10#${nsoill_out})), 'geogrid_file_input_grid': ${geogrid_file_input_grid}, - 'replace_vgtyp': ${replace_vgtyp}, - 'replace_sotyp': ${replace_sotyp}, - 'replace_vgfrc': ${replace_vgfrc}, + 'vgtyp_from_climo': ${vgtyp_from_climo}, + 'sotyp_from_climo': ${sotyp_from_climo}, + 'vgfrc_from_climo': ${vgfrc_from_climo}, + 'minmax_vgfrc_from_climo': ${minmax_vgfrc_from_climo}, + 'lai_from_climo': ${lai_from_climo}, 'tg3_from_soil': ${tg3_from_soil}, + 'thomp_mp_climo_file': ${thomp_mp_climo_file}, } " # @@ -633,7 +607,7 @@ mv_vrfy out.sfc.tile${TILE_RGNL}.nc \ mv_vrfy gfs_ctrl.nc ${ics_dir} -mv_vrfy gfs_bndy.nc ${ics_dir}/gfs_bndy.tile${TILE_RGNL}.000.nc +mv_vrfy gfs.bndy.nc ${ics_dir}/gfs_bndy.tile${TILE_RGNL}.000.nc # #----------------------------------------------------------------------- # diff --git a/scripts/exregional_make_lbcs.sh b/scripts/exregional_make_lbcs.sh index a284aa18b..84f317d74 100755 --- a/scripts/exregional_make_lbcs.sh +++ b/scripts/exregional_make_lbcs.sh @@ -57,7 +57,6 @@ hour zero). # valid_args=( \ "lbcs_dir" \ -"wgrib2_dir" \ "APRUN" \ ) process_args valid_args "$@" @@ -95,34 +94,31 @@ cd_vrfy $workdir # #----------------------------------------------------------------------- # -# Set physics-suite-dependent variables that are needed in the FORTRAN -# namelist file that the chgres executable will read in. +# Set physics-suite-dependent variable mapping table needed in the FORTRAN +# namelist file that the chgres_cube executable will read in. # #----------------------------------------------------------------------- # -phys_suite="" +varmap_file="" case "${CCPP_PHYS_SUITE}" in -"FV3_GFS_2017_gfdlmp" | "FV3_GFS_2017_gfdlmp_regional") - phys_suite="GFS" +"FV3_GFS_2017_gfdlmp" | "FV3_GFS_2017_gfdlmp_regional" | "FV3_GFS_v16beta" | \ +"FV3_GFS_v15p2" | "FV3_CPT_v0" ) + varmap_file="GFSphys_var_map.txt" ;; -"FV3_GSD_v0" | "FV3_GSD_SAR" | "FV3_RRFS_v1beta" ) - phys_suite="GSD" - ;; -"FV3_CPT_v0" ) - phys_suite="CPT" - ;; -"FV3_GFS_v15p2" ) - phys_suite="v15p2" - ;; -"FV3_GFS_v16beta" ) - phys_suite="v16beta" +"FV3_GSD_v0" | "FV3_GSD_SAR" | \ +"FV3_RRFS_v1beta" ) + if [ "${EXTRN_MDL_NAME_LBCS}" = "RAPX" ] || [ "${EXTRN_MDL_NAME_LBCS}" = "HRRRX" ]; then + varmap_file="GSDphys_var_map.txt" + elif [ "${EXTRN_MDL_NAME_LBCS}" = "NAM" ] || [ "${EXTRN_MDL_NAME_LBCS}" = "FV3GFS" ] || \ + [ "${EXTRN_MDL_NAME_LBCS}" = "GSMGFS" ]; then + varmap_file="GFSphys_var_map.txt" + fi ;; *) print_err_msg_exit "\ -Physics-suite-dependent namelist variables have not yet been specified -for this physics suite: +A variable mapping table has not yet been defined for this physics suite: CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" ;; @@ -131,7 +127,7 @@ esac #----------------------------------------------------------------------- # # Set external-model-dependent variables that are needed in the FORTRAN -# namelist file that the chgres executable will read in. These are de- +# namelist file that the chgres_cube executable will read in. These are de- # scribed below. Note that for a given external model, usually only a # subset of these all variables are set (since some may be irrelevant). # @@ -139,15 +135,19 @@ esac # Name of the external model from which we are obtaining the fields # needed to generate the LBCs. # +# fn_atm_nemsio: +# Name (not including path) of the nemsio file generated by the external +# model that contains the atmospheric fields. +# # fn_sfc_nemsio: # Name (not including path) of the nemsio file generated by the external # model that contains the surface fields. # # input_type: -# The "type" of input being provided to chgres. This contains a combi- +# The "type" of input being provided to chgres_cube. This contains a combi- # nation of information on the external model, external model file for- # mat, and maybe other parameters. For clarity, it would be best to -# eliminate this variable in chgres and replace with with 2 or 3 others +# eliminate this variable in chgres_cube and replace with with 2 or 3 others # (e.g. extrn_mdl, extrn_mdl_file_format, etc). # # tracers_input: @@ -163,30 +163,12 @@ esac # ment of tracers should be the name to use for the O3 mixing ratio in # the output file. For GSD physics, three additional tracers -- ice, # rain, and water number concentrations -- may be specified at the end -# of tracers, and these will be calculated by chgres. -# -# numsoil_out: -# The number of soil layers to include in the output NetCDF file. -# -# replace_FIELD, where FIELD="vgtyp", "sotyp", or "vgfrc": -# Logical variable indicating whether or not to obtain the field in -# question from climatology instead of the external model. The field in -# question is one of vegetation type (FIELD="vgtyp"), soil type (FIELD= -# "sotyp"), and vegetation fraction (FIELD="vgfrc"). If replace_FIELD -# is set to ".true.", then the field is obtained from climatology (re- -# gardless of whether or not it exists in an external model file). If -# it is set to ".false.", then the field is obtained from the external -# model. If the external model file does not provide this field, then -# chgres prints out an error message and stops. -# -# tg3_from_soil: -# Logical variable indicating whether or not to set the tg3 soil tempe- # Needs to be verified. -# rature field to the temperature of the deepest soil layer. +# of tracers, and these will be calculated by chgres_cube. # #----------------------------------------------------------------------- # -# GSK comments about chgres: +# GSK comments about chgres_cube: # # The following are the three atmsopheric tracers that are in the atmo- # spheric analysis (atmanl) nemsio file for CDATE=2017100700: @@ -235,13 +217,7 @@ fn_grib2="" input_type="" tracers_input="\"\"" tracers="\"\"" -numsoil_out="" -#geogrid_file_input_grid="\"\"" -replace_vgtyp="" -replace_sotyp="" -replace_vgfrc="" -tg3_from_soil="" - +thomp_mp_climo_file="" case "${EXTRN_MDL_NAME_LBCS}" in @@ -250,118 +226,108 @@ case "${EXTRN_MDL_NAME_LBCS}" in external_model="GSMGFS" - input_type="gfs_gaussian" # For spectral GFS Gaussian grid in nemsio format. + input_type="gfs_gaussian_nemsio" # For spectral GFS Gaussian grid in nemsio format. tracers_input="[\"spfh\",\"clwmr\",\"o3mr\"]" tracers="[\"sphum\",\"liq_wat\",\"o3mr\"]" +# +# Use Thompson climatology for ice- and water-friendly aerosols if CCPP suite uses Thompson MP +# + if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then + thomp_mp_climo_file="" + elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then + thomp_mp_climo_file="${FIXam}/Thompson_MP_MONTHLY_CLIMO.nc" + else + print_err_msg_exit "\ + The chosen CCPP physics suite is unsupported at this time: + CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" + fi - numsoil_out="4" - replace_vgtyp=".true." - replace_sotyp=".true." - replace_vgfrc=".true." - tg3_from_soil=".false." ;; "FV3GFS") +# +# Use Thompson climatology for ice- and water-friendly aerosols if CCPP suite uses Thompson MP +# + if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then + thomp_mp_climo_file="" + elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then + thomp_mp_climo_file="${FIXam}/Thompson_MP_MONTHLY_CLIMO.nc" + else + print_err_msg_exit "\ + The chosen CCPP physics suite is unsupported at this time: + CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" + fi + if [ "${FV3GFS_FILE_FMT_LBCS}" = "nemsio" ]; then external_model="FV3GFS" - input_type="gaussian" # For FV3-GFS Gaussian grid in nemsio format. + input_type="gaussian_nemsio" # For FV3-GFS Gaussian grid in nemsio format. tracers_input="[\"spfh\",\"clwmr\",\"o3mr\",\"icmr\",\"rwmr\",\"snmr\",\"grle\"]" -# -# If CCPP is being used, then the list of atmospheric tracers to include -# in the output file depends on the physics suite. Hopefully, this me- -# thod of specifying output tracers will be replaced with a variable -# table (which should be specific to each combination of external model, -# external model file type, and physics suite). -# - if [ "${USE_CCPP}" = "TRUE" ]; then - if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then - tracers="[\"sphum\",\"liq_wat\",\"o3mr\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\"]" - elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then -# For GSD physics, add three additional tracers (the ice, rain and water -# number concentrations) that are required for Thompson microphysics. - tracers="[\"sphum\",\"liq_wat\",\"o3mr\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\",\"ice_nc\",\"rain_nc\",\"water_nc\"]" - else - print_err_msg_exit "\ -The parameter \"tracers\" has not been defined for this combination of -external model (EXTRN_MDL_NAME_LBCS), physics suite (CCPP_PHYS_SUITE), -and FV3GFS file type (FV3GFS_FILE_FMT_LBCS): - EXTRN_MDL_NAME_LBCS = \"${EXTRN_MDL_NAME_LBCS}\" - CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\" - FV3GFS_FILE_FMT_LBCS = \"${FV3GFS_FILE_FMT_LBCS}\"" - fi -# -# If CCPP is not being used, the only physics suite that can be used is -# GFS. -# - else - tracers="[\"sphum\",\"liq_wat\",\"o3mr\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\"]" - fi + tracers="[\"sphum\",\"liq_wat\",\"o3mr\",\"ice_wat\",\"rainwat\",\"snowwat\",\"graupel\"]" elif [ "${FV3GFS_FILE_FMT_LBCS}" = "grib2" ]; then external_model="GFS" - input_type="grib2" fn_grib2="${EXTRN_MDL_FNS[0]}" + input_type="grib2" fi - numsoil_out="4" - replace_vgtyp=".true." - replace_sotyp=".true." - replace_vgfrc=".true." - tg3_from_soil=".false." - ;; "RAPX") external_model="RAP" - input_type="grib2" - if [ "${USE_CCPP}" = "TRUE" ]; then - if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ + ;; + + +"NAM") + + external_model="NAM" + input_type="grib2" +# +# Use Thompson climatology for ice- and water-friendly aerosols if CCPP suite uses Thompson MP +# + if [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_2017_gfdlmp_regional" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_CPT_v0" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ + [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then + thomp_mp_climo_file="" + elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ [ "${CCPP_PHYS_SUITE}" = "FV3_RRFS_v1beta" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v15p2" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GFS_v16beta" ]; then - numsoil_out="4" - elif [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_v0" ] || \ - [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then - numsoil_out="9" - else - print_err_msg_exit "\ -The parameter \"numsoil_out\" has not been defined for this combination -of external model (EXTRN_MDL_NAME_LBCS) and physics suite (CCPP_PHYS_SUITE): - EXTRN_MDL_NAME_LBCS = \"${EXTRN_MDL_NAME_LBCS}\" - CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" - fi + [ "${CCPP_PHYS_SUITE}" = "FV3_GSD_SAR" ]; then + thomp_mp_climo_file="${FIXam}/Thompson_MP_MONTHLY_CLIMO.nc" + else + print_err_msg_exit "\ + The chosen CCPP physics suite is unsupported at this time: + CCPP_PHYS_SUITE = \"${CCPP_PHYS_SUITE}\"" fi - replace_vgtyp=".false." - replace_sotyp=".false." - replace_vgfrc=".false." - tg3_from_soil=".true." - ;; - *) print_err_msg_exit "\ External-model-dependent namelist variables have not yet been specified @@ -390,7 +356,7 @@ fi # #----------------------------------------------------------------------- # -# Loop through the LBC update times and run chgres for each such time to +# Loop through the LBC update times and run chgres_cube for each such time to # obtain an LBC file for each that can be used as input to the FV3-LAM. # #----------------------------------------------------------------------- @@ -403,7 +369,7 @@ for (( i=0; i<${num_fhrs}; i++ )); do fhr="${EXTRN_MDL_LBC_SPEC_FHRS[$i]}" # # Set external model output file name and file type/format. Note that -# these are now inputs into chgres. +# these are now inputs into chgres_cube. # fn_atm_nemsio="" fn_grib2="" @@ -420,11 +386,17 @@ for (( i=0; i<${num_fhrs}; i++ )); do fi ;; "RAPX") + fn_grib2="${EXTRN_MDL_FNS[$i]}" + ;; + "HRRRX") + fn_grib2="${EXTRN_MDL_FNS[$i]}" + ;; + "NAM") fn_grib2="${EXTRN_MDL_FNS[$i]}" ;; *) print_err_msg_exit "\ -The external model output file name to use in the chgres FORTRAN name- +The external model output file name to use in the chgres_cube FORTRAN name- list file has not specified for this external model: EXTRN_MDL_NAME_LBCS = \"${EXTRN_MDL_NAME_LBCS}\"" ;; @@ -450,9 +422,6 @@ list file has not specified for this external model: # # Build the FORTRAN namelist file that chgres_cube will read in. # -# QUESTION: -# Do numsoil_out, ..., tg3_from_soil need to be in this namelist (as -# they are for the ICs namelist)? # # Create a multiline variable that consists of a yaml-compliant string @@ -461,35 +430,32 @@ list file has not specified for this external model: # this variable will be passed to a python script that will create the # namelist file. # - settings=" +settings=" 'config': { + 'fix_dir_input_grid': ${FIXgsm}, 'fix_dir_target_grid': ${FIXLAM}, - 'mosaic_file_target_grid': ${FIXLAM}/${CRES}${DOT_OR_USCORE}mosaic.halo${NH4}.nc, + 'mosaic_file_target_grid': ${FIXLAM}/${CRES}${DOT_OR_USCORE}mosaic.halo$((10#${NH4})).nc, 'orog_dir_target_grid': ${FIXLAM}, - 'orog_files_target_grid': ${CRES}${DOT_OR_USCORE}oro_data.tile7.halo${NH4}.nc, + 'orog_files_target_grid': ${CRES}${DOT_OR_USCORE}oro_data.tile${TILE_RGNL}.halo$((10#${NH4})).nc, 'vcoord_file_target_grid': ${FIXam}/global_hyblev.l65.txt, - 'mosaic_file_input_grid': '', - 'orog_dir_input_grid': '', - 'base_install_dir': ${CHGRES_DIR}, - 'wgrib2_path': ${wgrib2_dir}, + 'varmap_file': ${UFS_UTILS_DIR}/parm/varmap_tables/${varmap_file}, 'data_dir_input_grid': ${extrn_mdl_staging_dir}, 'atm_files_input_grid': ${fn_atm_nemsio}, 'sfc_files_input_grid': ${fn_sfc_nemsio}, 'grib2_file_input_grid': \"${fn_grib2}\", - 'cycle_mon': $((10#$mm)), - 'cycle_day': $((10#$dd)), - 'cycle_hour': $((10#$hh)), + 'cycle_mon': $((10#${mm})), + 'cycle_day': $((10#${dd})), + 'cycle_hour': $((10#${hh})), 'convert_atm': True, - 'convert_sfc': False, - 'convert_nst': False, 'regional': 2, - 'halo_bndy': ${NH4}, + 'halo_bndy': $((10#${NH4})), + 'halo_blend': $((10#${HALO_BLEND})), 'input_type': ${input_type}, 'external_model': ${external_model}, 'tracers_input': ${tracers_input}, 'tracers': ${tracers}, - 'phys_suite': ${phys_suite}, - } + 'thomp_mp_climo_file': ${thomp_mp_climo_file}, +} " # # Call the python script to create the namelist file. @@ -538,7 +504,7 @@ located in the following directory: # that of the external model since their start times may be offset). # fcst_hhh_FV3LAM=$( printf "%03d" "${LBC_SPEC_FCST_HRS[$i]}" ) - mv_vrfy gfs_bndy.nc ${lbcs_dir}/gfs_bndy.tile7.${fcst_hhh_FV3LAM}.nc + mv_vrfy gfs.bndy.nc ${lbcs_dir}/gfs_bndy.tile7.${fcst_hhh_FV3LAM}.nc done # diff --git a/sorc/global_equiv_resol.fd/Makefile_cheyenne b/sorc/global_equiv_resol.fd/Makefile_cheyenne deleted file mode 100644 index 456bb3fed..000000000 --- a/sorc/global_equiv_resol.fd/Makefile_cheyenne +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a ${NETCDF}/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -EXEC = global_equiv_resol - -.PHONY: all -all : $(EXEC) - -$(EXEC): global_equiv_resol.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/global_equiv_resol.fd/Makefile_hera b/sorc/global_equiv_resol.fd/Makefile_hera deleted file mode 100644 index c1f884090..000000000 --- a/sorc/global_equiv_resol.fd/Makefile_hera +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -EXEC = global_equiv_resol - -.PHONY: all -all : $(EXEC) - -$(EXEC): global_equiv_resol.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/global_equiv_resol.fd/Makefile_jet b/sorc/global_equiv_resol.fd/Makefile_jet deleted file mode 100644 index c1f884090..000000000 --- a/sorc/global_equiv_resol.fd/Makefile_jet +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -EXEC = global_equiv_resol - -.PHONY: all -all : $(EXEC) - -$(EXEC): global_equiv_resol.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/global_equiv_resol.fd/Makefile_odin b/sorc/global_equiv_resol.fd/Makefile_odin deleted file mode 100644 index 9a6073837..000000000 --- a/sorc/global_equiv_resol.fd/Makefile_odin +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a -lz - -FC = ftn -FFLAGS = -g -O2 $(INC) - -EXEC = global_equiv_resol - -.PHONY: all -all : $(EXEC) - -$(EXEC): global_equiv_resol.o - $(FC) $(FFLAGS) -o $@ $^ $(LIBS) - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/global_equiv_resol.fd/Makefile_stampede b/sorc/global_equiv_resol.fd/Makefile_stampede deleted file mode 100644 index ee9355051..000000000 --- a/sorc/global_equiv_resol.fd/Makefile_stampede +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a ${HDF5}/lib/libsz.a -lz - -FC = mpiifort -FFLAGS = -g -O2 $(INC) - -EXEC = global_equiv_resol - -.PHONY: all -all : $(EXEC) - -$(EXEC): global_equiv_resol.o - $(FC) $(FFLAGS) -o $@ $^ $(LIBS) - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/global_equiv_resol.fd/Makefile_theia b/sorc/global_equiv_resol.fd/Makefile_theia deleted file mode 100644 index c1f884090..000000000 --- a/sorc/global_equiv_resol.fd/Makefile_theia +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -EXEC = global_equiv_resol - -.PHONY: all -all : $(EXEC) - -$(EXEC): global_equiv_resol.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/global_equiv_resol.fd/compile.sh b/sorc/global_equiv_resol.fd/compile.sh deleted file mode 100755 index 7edd72c3b..000000000 --- a/sorc/global_equiv_resol.fd/compile.sh +++ /dev/null @@ -1,9 +0,0 @@ -#!/bin/bash - -module purge -module load intel/18.1.163 -module load netcdf/4.6.1 -module load hdf5/1.10.4 -module list - -make diff --git a/sorc/global_equiv_resol.fd/global_equiv_resol.f90 b/sorc/global_equiv_resol.fd/global_equiv_resol.f90 deleted file mode 100644 index 7bd9c813f..000000000 --- a/sorc/global_equiv_resol.fd/global_equiv_resol.f90 +++ /dev/null @@ -1,185 +0,0 @@ -!======================================================================= -program global_equiv_resol -!======================================================================= - - use netcdf - - implicit none - - integer, parameter :: dp = kind(1.0d0) - real(dp), parameter :: pi_geom = 4.0*atan(1.0), & - radius_Earth = 6371000.0 - - character(len=256) :: grid_fn - integer :: ncid, nxSG_dimid, nySG_dimid, dASG_varid, num_args - integer :: nxSG, nySG, nx, ny, RES_equiv - real(dp) :: avg_cell_size, min_cell_size, max_cell_size - real(dp), dimension(:,:), allocatable :: & - quarter_dA_ll, quarter_dA_lr, quarter_dA_ur, quarter_dA_ul, & - dASG, dA, sqrt_dA -! -!======================================================================= -! -! Read in the name of the file from the command line. The command-line -! call to this program should have exactly one argument consisting of -! the path to the NetCDF grid specification file to be read in. If this -! is not the case, print out a usage message and exit. -! -!======================================================================= -! - num_args = command_argument_count() - if (num_args == 1) then - call get_command_argument(1, grid_fn) - else - WRITE(*,500) - WRITE(*,500) "Exactly one argument must be specified to program global_equiv_resol." - WRITE(*,500) "Usage:" - WRITE(*,500) - WRITE(*,500) " global_equiv_resol path_to_grid_file" - WRITE(*,500) - WRITE(*,500) "where path_to_grid_file is the path to the NetCDF grid file. Actual " - WRITE(*,500) "number of specified command line arguments is:" - WRITE(*,510) " num_args = ", num_args - WRITE(*,500) "Stopping." -500 FORMAT(A) -510 FORMAT(A, I3) - STOP - end if -! -!======================================================================= -! -! Open the grid file and read in the dimensions of the supergrid. The -! supergrid is a grid that has twice the resolution of the actual/compu- -! tational grid. In the file, the names of the supergrid dimensions are -! nx and ny. Here, however, we reserve those names for the dimensions -! of the actual grid (since in the FV3 code and in other data files, nx -! and ny are used to denote the dimensions of the actual grid) and in- -! stead use the variables nxSG and nySG to denote the dimensions of the -! supergrid. -! -!======================================================================= -! - WRITE(*,500) - WRITE(*,500) "Opening NetCDF grid file for reading/writing:" - WRITE(*,500) " grid_fn = " // trim(grid_fn) - - call check( nf90_open(trim(grid_fn), NF90_WRITE, ncid) ) - - call check( nf90_inq_dimid(ncid, "nx", nxSG_dimid) ) - call check( nf90_inquire_dimension(ncid, nxSG_dimid, len=nxSG) ) - - call check( nf90_inq_dimid(ncid, "ny", nySG_dimid) ) - call check( nf90_inquire_dimension(ncid, nySG_dimid, len=nySG) ) - - WRITE(*,500) - WRITE(*,500) "Dimensions of supergrid are:" - WRITE(*,520) " nxSG = ", nxSG - WRITE(*,520) " nySG = ", nySG -520 FORMAT(A, I7) -! -!======================================================================= -! -! Read in the cell areas on the supergrid. Then add the areas of the -! four supergrid cells that make up one grid cell to obtain the cell -! areas on the actual grid. -! -!======================================================================= -! - allocate(dASG(0:nxSG-1, 0:nySG-1)) - call check( nf90_inq_varid(ncid, "area", dASG_varid) ) - call check( nf90_get_var(ncid, dASG_varid, dASG) ) - - nx = nxSG/2 - ny = nySG/2 - - WRITE(*,500) - WRITE(*,500) "Dimensions of (actual, i.e. computational) grid are:" - WRITE(*,520) " nx = ", nx - WRITE(*,520) " ny = ", ny - - allocate(quarter_dA_ll(0:nx-1, 0:ny-1)) - allocate(quarter_dA_lr(0:nx-1, 0:ny-1)) - allocate(quarter_dA_ul(0:nx-1, 0:ny-1)) - allocate(quarter_dA_ur(0:nx-1, 0:ny-1)) - - quarter_dA_ll = dASG(0:nxSG-1:2, 0:nySG-1:2) - quarter_dA_lr = dASG(0:nxSG-1:2, 1:nySG-1:2) - quarter_dA_ur = dASG(1:nxSG-1:2, 1:nySG-1:2) - quarter_dA_ul = dASG(1:nxSG-1:2, 0:nySG-1:2) - - allocate(dA(0:nx-1, 0:ny-1)) - allocate(sqrt_dA(0:nx-1, 0:ny-1)) - - dA = quarter_dA_ll + quarter_dA_lr + quarter_dA_ur + quarter_dA_ul -! -!======================================================================= -! -! Calculate a typical/representative cell size for each cell by taking -! the square root of the area of the cell. Then calculate the minimum, -! maximum, and average cell sizes over the whole grid. -! -!======================================================================= -! - sqrt_dA = sqrt(dA) - min_cell_size = minval(sqrt_dA) - max_cell_size = maxval(sqrt_dA) - avg_cell_size = sum(sqrt_dA)/(nx*ny) - - WRITE(*,500) - WRITE(*,500) "Minimum, maximum, and average cell sizes are (based on square" - WRITE(*,500) "root of cell area):" - WRITE(*,530) " min_cell_size = ", min_cell_size - WRITE(*,530) " max_cell_size = ", max_cell_size - WRITE(*,530) " avg_cell_size = ", avg_cell_size -530 FORMAT(A, G) -! -!======================================================================= -! -! Use the average cell size to calculate an equivalent global uniform -! cubed-sphere resolution (in units of number of cells) for the regional -! grid. This is the RES that a global uniform (i.e. stretch factor of -! 1) cubed-sphere grid would need to have in order to have the same no- -! minal cell size as the average cell size of the regional grid. -! -!======================================================================= -! - RES_equiv = nint( (2.0*pi_geom*radius_Earth)/(4.0*avg_cell_size) ) - - WRITE(*,500) - WRITE(*,500) "Equivalent global uniform cubed-sphere resolution is:" - WRITE(*,530) " RES_equiv = ", RES_equiv -! -!======================================================================= -! -! Write the average cell size and equivalent global resolution to the -! grid file as a global attributes. -! -!======================================================================= -! - WRITE(*,500) - WRITE(*,500) "Writing avg_cell_size and RES_equiv to the grid specification" - WRITE(*,500) "file as global attributes..." - - call check( nf90_redef(ncid) ) - call check( nf90_put_att(ncid, NF90_GLOBAL, "avg_cell_size", avg_cell_size) ) - call check( nf90_put_att(ncid, NF90_GLOBAL, "RES_equiv", RES_equiv) ) - call check( nf90_enddef(ncid) ) - - call check( nf90_close(ncid) ) - - WRITE(*,500) - WRITE(*,500) "Done." - -end program global_equiv_resol - - -subroutine check(status) - use netcdf - integer,intent(in) :: status -! - if(status /= nf90_noerr) then - write(0,*) ' check netcdf status = ', status - write(0,'("error ", a)') trim(nf90_strerror(status)) - stop "Stopped" - endif -end subroutine check diff --git a/sorc/mosaic_file.fd/Makefile_cheyenne b/sorc/mosaic_file.fd/Makefile_cheyenne deleted file mode 100644 index b7d132ffc..000000000 --- a/sorc/mosaic_file.fd/Makefile_cheyenne +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a ${NETCDF}/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -EXEC = mosaic_file - -.PHONY: all -all : $(EXEC) - -$(EXEC): mosaic_file.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/mosaic_file.fd/Makefile_hera b/sorc/mosaic_file.fd/Makefile_hera deleted file mode 100644 index 9a405799c..000000000 --- a/sorc/mosaic_file.fd/Makefile_hera +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -EXEC = mosaic_file - -.PHONY: all -all : $(EXEC) - -$(EXEC): mosaic_file.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/mosaic_file.fd/Makefile_jet b/sorc/mosaic_file.fd/Makefile_jet deleted file mode 100644 index 9a405799c..000000000 --- a/sorc/mosaic_file.fd/Makefile_jet +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -EXEC = mosaic_file - -.PHONY: all -all : $(EXEC) - -$(EXEC): mosaic_file.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/mosaic_file.fd/Makefile_odin b/sorc/mosaic_file.fd/Makefile_odin deleted file mode 100644 index 97c191a3b..000000000 --- a/sorc/mosaic_file.fd/Makefile_odin +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a -lz - -FC = ftn -FFLAGS = -g -O2 $(INC) - -EXEC = mosaic_file - -.PHONY: all -all : $(EXEC) - -$(EXEC): mosaic_file.o - $(FC) $(FFLAGS) -o $@ $^ $(LIBS) - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/mosaic_file.fd/Makefile_stampede b/sorc/mosaic_file.fd/Makefile_stampede deleted file mode 100644 index 8a3b1a83a..000000000 --- a/sorc/mosaic_file.fd/Makefile_stampede +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a ${HDF5}/lib/libsz.a -lz - -FC = mpiifort -FFLAGS = -g -O2 $(INC) - -EXEC = mosaic_file - -.PHONY: all -all : $(EXEC) - -$(EXEC): mosaic_file.o - $(FC) $(FFLAGS) -o $@ $^ $(LIBS) - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/mosaic_file.fd/Makefile_theia b/sorc/mosaic_file.fd/Makefile_theia deleted file mode 100644 index 9a405799c..000000000 --- a/sorc/mosaic_file.fd/Makefile_theia +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -EXEC = mosaic_file - -.PHONY: all -all : $(EXEC) - -$(EXEC): mosaic_file.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(EXEC) diff --git a/sorc/mosaic_file.fd/compile.sh b/sorc/mosaic_file.fd/compile.sh deleted file mode 100755 index 7edd72c3b..000000000 --- a/sorc/mosaic_file.fd/compile.sh +++ /dev/null @@ -1,9 +0,0 @@ -#!/bin/bash - -module purge -module load intel/18.1.163 -module load netcdf/4.6.1 -module load hdf5/1.10.4 -module list - -make diff --git a/sorc/mosaic_file.fd/mosaic_file.f90 b/sorc/mosaic_file.fd/mosaic_file.f90 deleted file mode 100644 index 7263e45bd..000000000 --- a/sorc/mosaic_file.fd/mosaic_file.f90 +++ /dev/null @@ -1,126 +0,0 @@ -!======================================================================= -program mosaic_file -!======================================================================= - - use netcdf - - implicit none - - integer, parameter:: dp=kind(1.0d0) - - integer, parameter :: ntiles = 1 - integer, parameter :: string = 255 - - logical :: file_exists - integer :: num_args, i, ncid, & - ntiles_dimid, string_dimid, & - mosaic_varid, gridlocation_varid, & - gridfiles_varid, gridtiles_varid - integer, dimension(1) :: dimids1D - integer, dimension(2) :: dimids2D - integer, dimension(ntiles) :: tile_inds - character(len=string) :: mosaic_fn, mosaic, gridlocation, CRES, tmp_str -! -!======================================================================= -! -! Create the grid mosaic file that FV3 expects to be present in the IN- -! PUT subdirectory of the run directory. -! -!======================================================================= -! - num_args = command_argument_count() - if (num_args == 1) then - call get_command_argument(1, CRES) - else - WRITE(*,500) - WRITE(*,500) "Exactly one argument must be specified to program mosaic_file." - WRITE(*,500) "Usage:" - WRITE(*,500) - WRITE(*,500) " mosaic_file CRES" - WRITE(*,500) - WRITE(*,500) "where CRES is the cubed-sphere grid resolution that will" - WRITE(*,500) "be used to form the name of the grid specification file(s)" - WRITE(*,500) "stored in the variable gridfiles in the grid mosaic file." - WRITE(*,500) "Actual number of specified command line arguments is:" - WRITE(*,510) " num_args = ", num_args - WRITE(*,500) "Stopping." -500 FORMAT(A) -510 FORMAT(A, I3) - STOP - end if -! -!======================================================================= -! -! Create the grid mosaic file that FV3 expects to be present in the IN- -! PUT subdirectory of the run directory. Then create dimensions, varia- -! bles, and attributes within it. -! -!======================================================================= -! - mosaic_fn = trim(CRES) // "_mosaic.nc" - - call check( nf90_create(mosaic_fn, NF90_64BIT_OFFSET, ncid) ) - call check( nf90_def_dim(ncid, "ntiles", ntiles, ntiles_dimid) ) - call check( nf90_def_dim(ncid, "string", string, string_dimid) ) - - dimids1D = (/ string_dimid /) - call check( nf90_def_var(ncid, "mosaic", NF90_CHAR, dimids1D, mosaic_varid) ) - call check( nf90_put_att(ncid, mosaic_varid, "standard_name", "grid_mosaic_spec") ) - call check( nf90_put_att(ncid, mosaic_varid, "children", "gridtiles")) - call check( nf90_put_att(ncid, mosaic_varid, "contact_regions", "contacts") ) - call check( nf90_put_att(ncid, mosaic_varid, "grid_descriptor", "") ) - - dimids1D = (/ string_dimid /) - call check( nf90_def_var(ncid, "gridlocation", NF90_CHAR, dimids1D, gridlocation_varid) ) - call check( nf90_put_att(ncid, gridlocation_varid, "standard_name", "grid_file_location") ) - - dimids2D = (/ string_dimid, ntiles_dimid /) - call check( nf90_def_var(ncid, "gridfiles", NF90_CHAR, dimids2D, gridfiles_varid) ) - call check( nf90_def_var(ncid, "gridtiles", NF90_CHAR, dimids2D, gridtiles_varid) ) - - call check( nf90_put_att(ncid, NF90_GLOBAL, "grid_version", "") ) - call check( nf90_put_att(ncid, NF90_GLOBAL, "code_version", "") ) - call check( nf90_put_att(ncid, NF90_GLOBAL, "history", "") ) - - call check( nf90_enddef(ncid) ) -! -!======================================================================= -! -! Assign values to variables in the grid mosaic file. The only one that -! seems to be read by the FV3 code (at least in regional mode) is the -! string variable "gridfiles" that contains the name of the grid speci- -! fication file. -! -!======================================================================= -! - mosaic = mosaic_fn - gridlocation = "/path/to/directory" - - call check( nf90_put_var(ncid, mosaic_varid, trim(mosaic)) ) - call check( nf90_put_var(ncid, gridlocation_varid, trim(gridlocation))) - - tile_inds(1) = 7 - do i=1, ntiles - write(tmp_str, 520) "tile", tile_inds(i) - call check( nf90_put_var(ncid, gridtiles_varid, trim(tmp_str), start=(/1,i/)) ) - tmp_str = trim(CRES) // "_grid." // trim(tmp_str) // ".nc" - call check( nf90_put_var(ncid, gridfiles_varid, trim(tmp_str), start=(/1,i/)) ) - end do -520 FORMAT(A, I1) -530 FORMAT(4A) - - call check( nf90_close(ncid) ) - -end program mosaic_file - - -subroutine check(status) - use netcdf - integer,intent(in) :: status -! - if(status /= nf90_noerr) then - write(0,*) ' check netcdf status = ', status - write(0,'("error ", a)') trim(nf90_strerror(status)) - stop "Stopped" - endif -end subroutine check diff --git a/sorc/regional_grid.fd/Makefile_cheyenne b/sorc/regional_grid.fd/Makefile_cheyenne deleted file mode 100644 index 95f68ca54..000000000 --- a/sorc/regional_grid.fd/Makefile_cheyenne +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a ${NETCDF}/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -REGIONAL_GRID = regional_grid - -.PHONY: all -all : $(REGIONAL_GRID) - -$(REGIONAL_GRID): pkind.o pietc.o pmat.o pmat4.o pmat5.o psym2.o gen_schmidt.o hgrid_ak.o regional_grid.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(REGIONAL_GRID) diff --git a/sorc/regional_grid.fd/Makefile_hera b/sorc/regional_grid.fd/Makefile_hera deleted file mode 100644 index cdcb21494..000000000 --- a/sorc/regional_grid.fd/Makefile_hera +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -REGIONAL_GRID = regional_grid - -.PHONY: all -all : $(REGIONAL_GRID) - -$(REGIONAL_GRID): pkind.o pietc.o pmat.o pmat4.o pmat5.o psym2.o gen_schmidt.o hgrid_ak.o regional_grid.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(REGIONAL_GRID) diff --git a/sorc/regional_grid.fd/Makefile_jet b/sorc/regional_grid.fd/Makefile_jet deleted file mode 100644 index cdcb21494..000000000 --- a/sorc/regional_grid.fd/Makefile_jet +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -REGIONAL_GRID = regional_grid - -.PHONY: all -all : $(REGIONAL_GRID) - -$(REGIONAL_GRID): pkind.o pietc.o pmat.o pmat4.o pmat5.o psym2.o gen_schmidt.o hgrid_ak.o regional_grid.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(REGIONAL_GRID) diff --git a/sorc/regional_grid.fd/Makefile_odin b/sorc/regional_grid.fd/Makefile_odin deleted file mode 100644 index 5b143a753..000000000 --- a/sorc/regional_grid.fd/Makefile_odin +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a -lz - -FC = ftn -FFLAGS = -g -O2 $(INC) - -REGIONAL_GRID = regional_grid - -.PHONY: all -all : $(REGIONAL_GRID) - -$(REGIONAL_GRID): pkind.o pietc.o pmat.o pmat4.o pmat5.o psym2.o gen_schmidt.o hgrid_ak.o regional_grid.o - $(FC) $(FFLAGS) -o $@ $^ $(LIBS) - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(REGIONAL_GRID) diff --git a/sorc/regional_grid.fd/Makefile_stampede b/sorc/regional_grid.fd/Makefile_stampede deleted file mode 100644 index 5f584611b..000000000 --- a/sorc/regional_grid.fd/Makefile_stampede +++ /dev/null @@ -1,29 +0,0 @@ - - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a ${HDF5}/lib/libsz.a -lz - -FC = mpiifort -FFLAGS = -g -O2 $(INC) - -REGIONAL_GRID = regional_grid - -.PHONY: all -all : $(REGIONAL_GRID) - -$(REGIONAL_GRID): pkind.o pietc.o pmat.o pmat4.o pmat5.o psym2.o gen_schmidt.o hgrid_ak.o regional_grid.o - $(FC) $(FFLAGS) -o $@ $^ $(LIBS) - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(REGIONAL_GRID) diff --git a/sorc/regional_grid.fd/Makefile_theia b/sorc/regional_grid.fd/Makefile_theia deleted file mode 100644 index cdcb21494..000000000 --- a/sorc/regional_grid.fd/Makefile_theia +++ /dev/null @@ -1,29 +0,0 @@ -SHELL := bash - -MAKEFLAGS += --warn-undefined-variables - -INC = -I${NETCDF}/include - -LIBS = ${NETCDF}/lib/libnetcdff.a ${NETCDF}/lib/libnetcdf.a \ - ${HDF5}/lib/libhdf5_hl.a ${HDF5}/lib/libhdf5.a /apps/szip/2.1/lib/libsz.a -lz - -FC = ifort -FFLAGS = -g -O2 $(INC) - -REGIONAL_GRID = regional_grid - -.PHONY: all -all : $(REGIONAL_GRID) - -$(REGIONAL_GRID): pkind.o pietc.o pmat.o pmat4.o pmat5.o psym2.o gen_schmidt.o hgrid_ak.o regional_grid.o $(LIBS) - $(FC) $(FFLAGS) -o $@ $^ - -.SUFFIXES: -.SUFFIXES: .f90 .o - -.f90.o: - $(FC) $(FFLAGS) -c $< - -.PHONY: clean -clean: - rm -f *.o *.mod $(REGIONAL_GRID) diff --git a/sorc/regional_grid.fd/compile.sh b/sorc/regional_grid.fd/compile.sh deleted file mode 100755 index f931fce87..000000000 --- a/sorc/regional_grid.fd/compile.sh +++ /dev/null @@ -1,9 +0,0 @@ -#!/bin/bash - -module purge -module load intel/18.1.163 -module load netcdf/4.6.1 -module load hdf5/1.10.4 -module list - -make -f Makefile_theia diff --git a/sorc/regional_grid.fd/gen_schmidt.f90 b/sorc/regional_grid.fd/gen_schmidt.f90 deleted file mode 100644 index 68fd2b6e2..000000000 --- a/sorc/regional_grid.fd/gen_schmidt.f90 +++ /dev/null @@ -1,477 +0,0 @@ -!============================================================================= -subroutine get_qqt(nxh,nyh,ncor,j0xy,p,q) -!============================================================================= -! Assume the grid to be mirror-symmetric across both medians, so that the -! computation of the quality diagnostic, Q, need only involve the positive -! quadrant of the grid. The norm associated with the definition of Q is the -! Frobenius norm (Q is the grid-mean of the squared-Frobenius norm of the -! log of the Gram matrix of the given distribution of jacobian matrices.) -!============================================================================= -use pkind, only: dp -use pietc, only: u0,u1,o2 -use pmat4, only: outer_product -use psym2 -implicit none -integer, intent(in ):: nxh,nyh,ncor -real(dp),dimension(3,2,0:nxh,0:nyh),intent(in ):: j0xy -real(dp),dimension(2,2), intent(inout):: p -real(dp), intent( out):: q -!----------------------------------------------------------------------------- -integer,parameter :: nit=5 -real(dp),parameter :: acrit=1.e-8,dpx=.0099 -real(dp),dimension(0:nxh,0:nyh) :: wxy -real(dp),dimension(3,2) :: j0,j -real(dp),dimension(2,2) :: el,pf,elp,elmean,g,ppx,pmx,ppy,pmy -real(dp),dimension(2) :: hess,grad -real(dp) :: anorm,q00,qpx,qmx,qpy,qmy,c,w -integer :: ix,iy,lx,ly,it -!============================================================================= -call get_wxy(nxh,nyh,ncor,wxy)! <- get 2D extended trapezoidal averaging wts -if(p(1,1)==u0)then; p=0; p(1,1)=u1; p(2,2)=u1; endif -! Iteratively calibrate preconditioner, p, to make elmean vanish: -anorm=1 -do it=1,nit - elmean=0 - q=0 - do iy=0,nyh; do ix=0,nxh - j0=j0xy(:,:,ix,iy); w=wxy(ix,iy) -! Precondition the Jacobian using latest iteration of P: - j=matmul(j0,p) -! Find the Gram matrix, G, implied by the column vectors of the new J: - g=matmul(transpose(j),j) -! Find the matrix logarithm, L = log(G), contrinutions to elmean and q: - call logsym2(g,el); el=el/2; elmean=elmean+w*el; q=q+w*sum(el**2) - enddo ; enddo - if(anormnit)then - print'("WARNING: In get_qqt, apparent failure of iteration to converge")' - read(*,*) -endif - -q00=q -ppx=p; ppx(1,1)=ppx(1,1)*(1+dpx);qpx=0 -pmx=p; pmx(1,1)=pmx(1,1)*(1-dpx);qmx=0 -ppy=p; ppy(2,2)=ppy(2,2)*(1+dpx);qpy=0 -pmy=p; pmy(2,2)=pmy(2,2)*(1-dpx);qmy=0 -do iy=0,nyh; do ix=0,nxh - j0=j0xy(:,:,ix,iy); w=wxy(ix,iy) - j=matmul(j0,ppx); g=matmul(transpose(j),j) - call logsym2(g,el); el=el/2; qpx=qpx+w*sum(el**2) - j=matmul(j0,pmx); g=matmul(transpose(j),j) - call logsym2(g,el); el=el/2; qmx=qmx+w*sum(el**2) - j=matmul(j0,ppy); g=matmul(transpose(j),j) - call logsym2(g,el); el=el/2; qpy=qpy+w*sum(el**2) - j=matmul(j0,pmy); g=matmul(transpose(j),j) - call logsym2(g,el); el=el/2; qmy=qmy+w*sum(el**2) -enddo; enddo -! Estimate a (diagonal) Hessian matrix and a gradient vector: -hess=(/ (qpx-2*q00+qmx)/dpx**2, (qpy-2*q00+qmy)/dpx**2 /) -grad=(/ (qpx-qmx)/(2*dpx) , (qpy-qmy)/(2*dpx) /) - -!!print'('' hessian components:'',t30,2(1x,e20.14))',hess !!!!!!! -!!print'('' grad components:'',t30,2(1x,e20.14))',grad !!!!!!! -! If the hessian is positive, polish the final p with a final Newton iteration: -if(hess(1)>0 .and. hess(2)>0.)then - c=u1-grad(1)/hess(1); p(:,1)=p(:,1)*c - c=u1-grad(2)/hess(2); p(:,2)=p(:,2)*c -endif - -! and calculate the new q. Keep it only if is numerically smaller than before: -q00=0 -do iy=0,nyh; do ix=0,nxh - j0=j0xy(:,:,ix,iy); w=wxy(ix,iy) - j=matmul(j0,p); g=matmul(transpose(j),j) - call logsym2(g,el); el=el/2; q00=q00+w*sum(el**2) -enddo; enddo -!!print'('' adjusted final q: '',e20.14)',q00 -if(q00nit)then - print'("WARNING: In get_qqt, apparent failure of iteration to converge")' - read(*,*) -endif - -q00=q -ppx=p; ppx(1,1)=ppx(1,1)*(1+dpx);qpx=0 -pmx=p; pmx(1,1)=pmx(1,1)*(1-dpx);qmx=0 -ppy=p; ppy(2,2)=ppy(2,2)*(1+dpx);qpy=0 -pmy=p; pmy(2,2)=pmy(2,2)*(1-dpx);qmy=0 -do iy=0,nyh; do ix=0,nxh - j0=j0xy(:,:,ix,iy); w=wxy(ix,iy) - j=matmul(j0,ppx); g=matmul(transpose(j),j) - call logsym2(g,el);el=el/2;qpx=qpx+w*(twc*sum(el**2)+tw*(el(1,1)+el(2,2))**2) - j=matmul(j0,pmx); g=matmul(transpose(j),j) - call logsym2(g,el);el=el/2;qmx=qmx+w*(twc*sum(el**2)+tw*(el(1,1)+el(2,2))**2) - j=matmul(j0,ppy); g=matmul(transpose(j),j) - call logsym2(g,el);el=el/2;qpy=qpy+w*(twc*sum(el**2)+tw*(el(1,1)+el(2,2))**2) - j=matmul(j0,pmy); g=matmul(transpose(j),j) - call logsym2(g,el);el=el/2;qmy=qmy+w*(twc*sum(el**2)+tw*(el(1,1)+el(2,2))**2) -enddo; enddo -! Estimate a (diagonal) Hessian matrix and a gradient vector: -hess=(/ (qpx-2*q00+qmx)/dpx**2, (qpy-2*q00+qmy)/dpx**2 /) -hess=(/8.,8./) -grad=(/ (qpx-qmx)/(2*dpx) , (qpy-qmy)/(2*dpx) /) -! If the hessian is positive, polish p with a final Newton iteration: -if(hess(1)>0 .and. hess(2)>0.)then - c=u1-grad(1)/hess(1); p(:,1)=p(:,1)*c - c=u1-grad(2)/hess(2); p(:,2)=p(:,2)*c -endif - -! and calculate the new q. Keep it only if it's numerically smaller than before: -q00=0 -do iy=0,nyh; do ix=0,nxh - j0=j0xy(:,:,ix,iy); w=wxy(ix,iy) - j=matmul(j0,p); g=matmul(transpose(j),j) - call logsym2(g,el);el=el/2;q00=q00+w*(twc*sum(el**2)+tw*(el(1,1)+el(2,2))**2) -enddo; enddo -if(q004)stop 'In get_wxy; ncor is out of bounds' -if(ncor>=min(nxh,nyh))stop 'In get_wxy; ncor is too large for this small grid' -! the wx and wy are the weight coefficients for an unnormalized -! extended trapezoidal integration. The end correction coefficients can -! be found by staggering, then summing, the Adams-Moulton coefficients -! at both ends. -wx=u1; wx(0)=o2; wx(nxh:nxh-ncor:-1)=cor -wy=u1; wy(0)=o2; wy(nyh:nyh-ncor:-1)=cor -wxy=outer_product(wx,wy); wxy=wxy/sum(wxy) -end subroutine get_wxy - -!============================================================================= -subroutine getedges(arcx,arcy,edgex,edgey) -!============================================================================= -! For angles (degrees) of the arcs spanning the halfwidths between the -! region's center and its x and y edges, get the two cartesian vectors -! that represent the locations of these edge midpoints in the positive x and y -! directions. -!============================================================================= -use pkind, only: dp -use pietc, only: u0,dtor -implicit none -real(dp), intent(in ):: arcx,arcy -real(dp),dimension(3),intent(out):: edgex,edgey -!------------------------------------------------------------------------------ -real(dp):: cx,sx,cy,sy -!============================================================================== -cx=cos(arcx*dtor); sx=sin(arcx*dtor) -cy=cos(arcy*dtor); sy=sin(arcy*dtor) -edgex=(/sx,u0,cx/); edgey=(/u0,sy,cy/) -end subroutine getedges - -!============================================================================== -subroutine xmtoxc_ak(a,kappa,xm,xc,xcd,ff) -!============================================================================== -! Assuming the A-kappa parameterization of the generalized schmidt-transformed -! gnomonic mapping, and given a map-space 2-vector, xm, find the corresponding -! cartesian unit 3-vector and its derivative wrt xm, jacobian matrix, xcd. -! If for any reason the mapping cannot be done, return a raised failure -! flag, FF. -!============================================================================= -use pkind, only: dp -use pietc, only: T,F -implicit none -real(dp), intent(in ):: a,kappa -real(dp),dimension(2), intent(in ):: xm -real(dp),dimension(3), intent(out):: xc -real(dp),dimension(3,2),intent(out):: xcd -logical, intent(out):: ff -!----------------------------------------------------------------------------- -real(dp),dimension(2,2):: xtd,xsd -real(dp),dimension(2) :: xt,xs -!============================================================================= -call xmtoxt(a,xm,xt,xtd,ff); if(ff)return -call xttoxs(kappa,xt,xs,xsd,ff); if(ff)return -xsd=matmul(xsd,xtd) -call xstoxc(xs,xc,xcd) -xcd=matmul(xcd,xsd) -end subroutine xmtoxc_ak - -!============================================================================= -subroutine xctoxm_ak(a,kappa,xc,xm,ff) -!============================================================================= -! Inverse mapping of xmtoxc_ak. That is, go from given cartesian unit 3-vector, -! xc, to map coordinate 2-vector xm (or return a raised failure flag, FF, if -! the attempt fails). -!============================================================================= -use pkind, only: dp -use pietc, only: F,T,u0,u1 -implicit none -real(dp), intent(in ):: a,kappa -real(dp),dimension(3),intent(in ):: xc -real(dp),dimension(2),intent(out):: xm -logical, intent(out):: ff -!----------------------------------------------------------------------------- -real(dp),dimension(2):: xs,xt -!============================================================================= -ff=F -call xctoxs(xc,xs) -call xstoxt(kappa,xs,xt,ff); if(ff)return -call xttoxm(a,xt,xm,ff) -end subroutine xctoxm_ak - -!============================================================================== -subroutine zmtozt(a,zm,zt,ztd,ff) -!============================================================================== -! Evaluate the function, zt = tan(sqrt(A)*z)/sqrt(A), and its deivative, ztd, -! for positive and negative A and for the limiting case, A --> 0 -!============================================================================== -use pkind, only: dp -use pietc, only: F,T,u1,pih -implicit none -real(dp),intent(in ):: a,zm -real(dp),intent(out):: zt,ztd -logical, intent(out):: ff -!------------------------------------------------------------------------------ -real(dp):: ra -!============================================================================== -ff=f -if (a>0)then; ra=sqrt( a); zt=tan (ra*zm)/ra; ff=abs(ra*zm)>=pih -elseif(a<0)then; ra=sqrt(-a); zt=tanh(ra*zm)/ra -else ; zt=zm -endif -ztd=u1+a*zt*zt -end subroutine zmtozt - -!============================================================================= -subroutine zttozm(a,zt,zm,ff) -!============================================================================= -! Inverse of zmtozt -!============================================================================= -use pkind, only: dp -use pietc, only: F,u1 -implicit none -real(dp),intent(in ):: a,zt -real(dp),intent(out):: zm -logical, intent(out):: ff -!----------------------------------------------------------------------------- -real(dp):: ra,razt -!============================================================================= -ff=F -if (a>0)then; ra=sqrt( a); razt=ra*zt; zm=atan (razt)/ra -elseif(a<0)then; ra=sqrt(-a); razt=ra*zt; ff=abs(razt)>=u1; if(ff)return - zm=atanh(razt)/ra -else ; zm=zt -endif -end subroutine zttozm - -!============================================================================== -subroutine xmtoxt(a,xm,xt,xtd,ff) -!============================================================================== -! Like zmtozt, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd -!============================================================================== -use pkind, only: dp -implicit none -real(dp), intent(in ):: a -real(dp),dimension(2), intent(in ):: xm -real(dp),dimension(2), intent(out):: xt -real(dp),dimension(2,2),intent(out):: xtd -logical, intent(out):: ff -!----------------------------------------------------------------------------- -integer:: i -!============================================================================== -xtd=0; do i=1,2; call zmtozt(a,xm(i),xt(i),xtd(i,i),ff); if(ff)return; enddo -end subroutine xmtoxt - -!============================================================================= -subroutine xttoxm(a,xt,xm,ff) -!============================================================================= -! Inverse of xmtoxt -!============================================================================ -use pkind, only: dp -use pietc, only: F -implicit none -real(dp), intent(in ):: a -real(dp),dimension(2),intent(in ):: xt -real(dp),dimension(2),intent(out):: xm -logical ,intent(out):: ff -!----------------------------------------------------------------------------- -integer:: i -!============================================================================= -do i=1,2; call zttozm(a,xt(i),xm(i),ff); if(ff)return; enddo -end subroutine xttoxm - -!============================================================================== -subroutine xttoxs(kappa,xt,xs,xsd,ff) -!============================================================================== -! Scaled gnomonic plane xt to standard stereographic plane xs -!============================================================================== -use pkind, only: dp -use pietc, only: u0,u1 -implicit none -real(dp), intent(in ):: kappa -real(dp),dimension(2), intent(in ):: xt -real(dp),dimension(2), intent(out):: xs -real(dp),dimension(2,2),intent(out):: xsd -logical, intent(out):: ff -!------------------------------------------------------------------------------ -real(dp):: s,sp,rsp,rspp,rspps,rspdx,rspdy -!============================================================================== -s=kappa*(xt(1)*xt(1) + xt(2)*xt(2)); sp=u1+s -ff=(sp<=u0); if(ff)return -rsp=sqrt(sp) -rspp=u1+rsp -rspps=rspp**2 -xs=xt/rspp -rspdx=kappa*xt(1)/rsp -rspdy=kappa*xt(2)/rsp -xsd(1,1)=u1/rspp -xt(1)*rspdx/rspps -xsd(1,2)= -xt(1)*rspdy/rspps -xsd(2,1)= -xt(2)*rspdx/rspps -xsd(2,2)=u1/rspp -xt(2)*rspdy/rspps -end subroutine xttoxs - -!============================================================================= -subroutine xstoxt(kappa,xs,xt,ff) -!============================================================================= -! Inverse of xttoxs. -!============================================================================= -use pkind, only: dp -use pietc, only: u1 -implicit none -real(dp), intent(in ):: kappa -real(dp),dimension(2),intent(in ):: xs -real(dp),dimension(2),intent(out):: xt -logical, intent(out):: ff -!----------------------------------------------------------------------------- -real(dp):: s,sc -!============================================================================= -s=kappa*(xs(1)*xs(1)+xs(2)*xs(2)); sc=u1-s -ff=(sc<=0); if(ff)return -xt=2*xs/sc -end subroutine xstoxt - -!============================================================================= -subroutine xstoxc(xs,xc,xcd) -!============================================================================= -! Standard transformation from polar stereographic map coordinates, xs, to -! cartesian, xc, assuming the projection axis is polar. -! xcd=d(xc)/d(xs) is the jacobian matrix -!============================================================================= -use pkind, only: dp -use pietc, only: u1,u2 -use pmat4, only: outer_product -implicit none -real(dp),dimension(2), intent(in ):: xs -real(dp),dimension(3), intent(out):: xc -real(dp),dimension(3,2),intent(out):: xcd -!----------------------------------------------------------------------------- -real(dp):: zp -!============================================================================= -zp=u2/(u1+dot_product(xs,xs)); xc(1:2)=xs*zp; xc(3)=zp -xcd=-outer_product(xc,xs)*zp; xcd(1,1)=xcd(1,1)+zp; xcd(2,2)=xcd(2,2)+zp -xc(3)=xc(3)-u1 -end subroutine xstoxc - -!============================================================================= -subroutine xctoxs(xc,xs) -!============================================================================= -! Inverse of xstoxc. I.e., cartesians to stereographic -!============================================================================= -use pkind, only: dp -use pietc, only: u1 -implicit none -real(dp),dimension(3),intent(in ):: xc -real(dp),dimension(2),intent(out):: xs -!----------------------------------------------------------------------------- -real(dp):: zp -!============================================================================= -zp=u1+xc(3); xs=xc(1:2)/zp -end subroutine xctoxs diff --git a/sorc/regional_grid.fd/hgrid_ak.f90 b/sorc/regional_grid.fd/hgrid_ak.f90 deleted file mode 100644 index 389f53a69..000000000 --- a/sorc/regional_grid.fd/hgrid_ak.f90 +++ /dev/null @@ -1,104 +0,0 @@ -!============================================================================= -subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & - re,delx,dely, glat,glon,garea, ff) -!============================================================================= -! Use a and k as the parameters of a generalized Schmidt-transformed -! gnomonic mapping centered at (plat,plon) and twisted about this center -! by an azimuth angle of pazi counterclockwise (these angles in radians). -! -! Assuming the radius of the earth is re, and using the central mapping -! point as the coordinate origin, set up the grid with central x-spacing delx -! and y-spacing dely in physical units, and with the location of the left-lower -! corner of the grid at center-relative grid index pair, (lx,ly) and with -! the number of the grid spaces in x and y directions given by nx and ny. -! (Note that, for a centered rectangular grid lx and ly are negative and, in -! magnitude, half the values of nx and ny respectively.) -! Return the latitude and longitude, in radians again, of the grid points thus -! defined in the arrays, glat and glon, and return a rectangular array, garea, -! of dimensions nx-1 by ny-1, that contains the areas of each of the grid cells -! in the SQUARE of the same physical length unit that was employed to define -! the radius of the earth, re (and the central grid steps, delx and dely). -! -! if all goes well, return a .FALSE. failure flag, ff. If, for some -! reason, it is not possible to complete this task, return the failure flag -! as .TRUE. -!============================================================================= -use pkind, only: dp -use pietc, only: u0,u1,dtor -use pmat4, only: sarea -use pmat5, only: ctog -implicit none -integer, intent(in ):: lx,ly,nx,ny -real(dp), intent(in ):: a,k,plat,plon,pazi, & - re,delx,dely -real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon -real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea -logical, intent(out):: ff -!----------------------------------------------------------------------------- -real(dp),dimension(3,3):: prot,azirot -real(dp),dimension(3,2):: xcd -real(dp),dimension(3) :: xc -real(dp),dimension(2) :: xm -real(dp) :: clat,slat,clon,slon,cazi,sazi,& - rlat,drlata,drlatb,drlatc, & - rlon,drlona,drlonb,drlonc, rre -integer :: ix,iy,mx,my -!============================================================================= -clat=cos(plat); slat=sin(plat) -clon=cos(plon); slon=sin(plon) -cazi=cos(pazi); sazi=sin(pazi) - -azirot(:,1)=(/ cazi, sazi, u0/) -azirot(:,2)=(/-sazi, cazi, u0/) -azirot(:,3)=(/ u0, u0, u1/) - -prot(:,1)=(/ -slon, clon, u0/) -prot(:,2)=(/-slat*clon, -slat*slon, clat/) -prot(:,3)=(/ clat*clon, clat*slon, slat/) -prot=matmul(prot,azirot) -mx=lx+nx ! Index of the 'right' edge of the rectangular grid -my=ly+ny ! Index of the 'top' edge of the rectangular grid -!This code assumes symmetry about the grid center -do iy=ly,my - !xm(2)=iy*dely/re - xm(2)=-iy*dely/re - do ix=lx,mx - !xm(1)=ix*delx/re - xm(1)=-ix*delx/re - call xmtoxc_ak(a,k,xm,xc,xcd,ff) - if(ff)return - xcd=matmul(prot,xcd) - xc =matmul(prot,xc ) - call ctog(xc,glat(ix,iy),glon(ix,iy)) - enddo -enddo - -! Convert degrees to radians in the glat and glon arrays: -glat=glat*dtor -glon=glon*dtor - -! Compute the areas of the quadrilateral grid cells: -do iy=ly,my-1 - do ix=lx,mx-1 - rlat =glat(ix ,iy ) - drlata=glat(ix+1,iy )-rlat - drlatb=glat(ix+1,iy+1)-rlat - drlatc=glat(ix ,iy+1)-rlat - rlon =glon(ix ,iy ) - drlona=glon(ix+1,iy )-rlon - drlonb=glon(ix+1,iy+1)-rlon - drlonc=glon(ix ,iy+1)-rlon -! If 'I' is the grid point (ix,iy), 'A' is (ix+1,iy); 'B' is (ix+1,iy+1) -! and 'C' is (ix,iy+1), then the area of the grid cell IABC is the sum of -! the areas of the traingles, IAB and IBC (the latter being the negative -! of the signed area of triangle, ICB): - garea(ix,iy)=sarea(rlat, drlata,drlona, drlatb,drlonb) & - -sarea(rlat, drlatc,drlonc, drlatb,drlonb) - enddo -enddo -! Convert the areas to area units consistent with the length units used for -! the radius, re, of the sphere: -rre=re*re -garea=garea*rre - -end subroutine hgrid_ak diff --git a/sorc/regional_grid.fd/pietc.f90 b/sorc/regional_grid.fd/pietc.f90 deleted file mode 100644 index f6abe1af3..000000000 --- a/sorc/regional_grid.fd/pietc.f90 +++ /dev/null @@ -1,95 +0,0 @@ -! -!============================================================================= -module pietc -!============================================================================= -! R. J. Purser (jim.purser@noaa.gov) 2014 -! Some of the commonly used constants (pi etc) mainly for double-precision -! subroutines. -! ms10 etc are needed to satisfy the some (eg., gnu fortran) compilers' -! more rigorous standards regarding the way "data" statements are initialized. -! Zero and the first few units are u0,u1,u2, etc., their reciprocals being, -! o2,o3 etc and their square roots, r2,r3. Reciprocal roots are or2,or3 etc. -!============================================================================= -use pkind, only: dp,dpc -implicit none -logical ,parameter:: T=.true.,F=.false. !<- for pain-relief in logical ops -real(dp),parameter:: & - u0=0,u1=1,mu1=-u1,u2=2,mu2=-u2,u3=3,mu3=-u3,u4=4,mu4=-u4,u5=5,mu5=-u5, & - u6=6,mu6=-u6,o2=u1/2,o3=u1/3,o4=u1/4,o5=u1/5,o6=u1/6, & - pi =3.1415926535897932384626433832795028841971693993751058209749e0_dp, & - pi2=6.2831853071795864769252867665590057683943387987502116419498e0_dp, & - pih=1.5707963267948966192313216916397514420985846996875529104874e0_dp, & - rpi=1.7724538509055160272981674833411451827975494561223871282138e0_dp, & -! Important square-roots - r2 =1.4142135623730950488016887242096980785696718753769480731766e0_dp, & - r3 =1.7320508075688772935274463415058723669428052538103806280558e0_dp, & - r5 =2.2360679774997896964091736687312762354406183596115257242708e0_dp, & - or2=u1/r2,or3=u1/r3,or5=u1/r5, & -! Golden number: - phi=1.6180339887498948482045868343656381177203091798057628621354e0_dp, & -! Euler-Mascheroni constant: - euler=0.57721566490153286060651209008240243104215933593992359880e0_dp, & -! Degree to radians; radians to degrees: - dtor=pi/180,rtod=180/pi, & -! Sines of all main fractions of 90 degrees (down to ninths): - s10=.173648177666930348851716626769314796000375677184069387236241e0_dp,& - s11=.195090322016128267848284868477022240927691617751954807754502e0_dp,& - s13=.222520933956314404288902564496794759466355568764544955311987e0_dp,& - s15=.258819045102520762348898837624048328349068901319930513814003e0_dp,& - s18=.309016994374947424102293417182819058860154589902881431067724e0_dp,& - s20=.342020143325668733044099614682259580763083367514160628465048e0_dp,& - s22=.382683432365089771728459984030398866761344562485627041433800e0_dp,& - s26=.433883739117558120475768332848358754609990727787459876444547e0_dp,& - s30=o2, & - s34=.555570233019602224742830813948532874374937190754804045924153e0_dp,& - s36=.587785252292473129168705954639072768597652437643145991072272e0_dp,& - s39=.623489801858733530525004884004239810632274730896402105365549e0_dp,& - s40=.642787609686539326322643409907263432907559884205681790324977e0_dp,& - s45=or2, & - s50=.766044443118978035202392650555416673935832457080395245854045e0_dp,& - s51=.781831482468029808708444526674057750232334518708687528980634e0_dp,& - s54=.809016994374947424102293417182819058860154589902881431067724e0_dp,& - s56=.831469612302545237078788377617905756738560811987249963446124e0_dp,& - s60=r3*o2, & - s64=.900968867902419126236102319507445051165919162131857150053562e0_dp,& - s68=.923879532511286756128183189396788286822416625863642486115097e0_dp,& - s70=.939692620785908384054109277324731469936208134264464633090286e0_dp,& - s72=.951056516295153572116439333379382143405698634125750222447305e0_dp,& - s75=.965925826289068286749743199728897367633904839008404550402343e0_dp,& - s77=.974927912181823607018131682993931217232785800619997437648079e0_dp,& - s79=.980785280403230449126182236134239036973933730893336095002916e0_dp,& - s80=.984807753012208059366743024589523013670643251719842418790025e0_dp,& -! ... and their minuses: - ms10=-s10,ms11=-s11,ms13=-s13,ms15=-s15,ms18=-s18,ms20=-s20,ms22=-s22,& - ms26=-s26,ms30=-s30,ms34=-s34,ms36=-s36,ms39=-s39,ms40=-s40,ms45=-s45,& - ms50=-s50,ms51=-s51,ms54=-s54,ms56=-s56,ms60=-s60,ms64=-s64,ms68=-s68,& - ms70=-s70,ms72=-s72,ms75=-s75,ms77=-s77,ms79=-s79,ms80=-s80 - -complex(dpc),parameter:: & - c0=(u0,u0),c1=(u1,u0),mc1=-c1,ci=(u0,u1),mci=-ci,cipi=ci*pi, & -! Main fractional rotations, as unimodular complex numbers: - z000=c1 ,z010=( s80,s10),z011=( s79,s11),z013=( s77,s13),& - z015=( s75,s15),z018=( s72,s18),z020=( s70,s20),z022=( s68,s22),& - z026=( s64,s26),z030=( s60,s30),z034=( s56,s34),z036=( s54,s36),& - z039=( s51,s39),z040=( s50,s40),z045=( s45,s45),z050=( s40,s50),& - z051=( s39,s51),z054=( s36,s54),z056=( s34,s56),z060=( s30,s60),& - z064=( s26,s64),z068=( s22,s68),z070=( s20,s70),z072=( s18,s72),& - z075=( s15,s75),z077=( s13,s77),z079=( s11,s79),z080=( s10,s80),& - z090=ci, z100=(ms10,s80),z101=(ms11,s79),z103=(ms13,s77),& - z105=(ms15,s75),z108=(ms18,s72),z110=(ms20,s70),z112=(ms22,s68),& - z116=(ms26,s64),z120=(ms30,s60),z124=(ms34,s56),z126=(ms36,s54),& - z129=(ms39,s51),z130=(ms40,s50),z135=(ms45,s45),z140=(ms50,s40),& - z141=(ms51,s39),z144=(ms54,s36),z146=(ms56,s34),z150=(ms60,s30),& - z154=(ms64,s26),z158=(ms68,s22),z160=(ms70,s20),z162=(ms72,s18),& - z165=(ms75,s15),z167=(ms77,s13),z169=(ms79,s11),z170=(ms80,s10),& - z180=-z000,z190=-z010,z191=-z011,z193=-z013,z195=-z015,z198=-z018,& - z200=-z020,z202=-z022,z206=-z026,z210=-z030,z214=-z034,z216=-z036,& - z219=-z039,z220=-z040,z225=-z045,z230=-z050,z231=-z051,z234=-z054,& - z236=-z056,z240=-z060,z244=-z064,z248=-z068,z250=-z070,z252=-z072,& - z255=-z075,z257=-z077,z259=-z079,z260=-z080,z270=-z090,z280=-z100,& - z281=-z101,z283=-z103,z285=-z105,z288=-z108,z290=-z110,z292=-z112,& - z296=-z116,z300=-z120,z304=-z124,z306=-z126,z309=-z129,z310=-z130,& - z315=-z135,z320=-z140,z321=-z141,z324=-z144,z326=-z146,z330=-z150,& - z334=-z154,z338=-z158,z340=-z160,z342=-z162,z345=-z165,z347=-z167,& - z349=-z169,z350=-z170 -end module pietc diff --git a/sorc/regional_grid.fd/pkind.f90 b/sorc/regional_grid.fd/pkind.f90 deleted file mode 100644 index abc5841b7..000000000 --- a/sorc/regional_grid.fd/pkind.f90 +++ /dev/null @@ -1,8 +0,0 @@ -module pkind -private:: one_dpi; integer(8),parameter:: one_dpi=1 -integer,parameter:: dpi=kind(one_dpi) -integer,parameter:: sp=kind(1.0) -integer,parameter:: dp=kind(1.0d0) -integer,parameter:: spc=kind((1.0,1.0)) -integer,parameter:: dpc=kind((1.0d0,1.0d0)) -end module pkind diff --git a/sorc/regional_grid.fd/pmat.f90 b/sorc/regional_grid.fd/pmat.f90 deleted file mode 100644 index 11f1b0f7b..000000000 --- a/sorc/regional_grid.fd/pmat.f90 +++ /dev/null @@ -1,1082 +0,0 @@ -! -! ********************************************** -! * MODULE pmat * -! * R. J. Purser, NOAA/NCEP/EMC 1993 * -! * and Tsukasa Fujita, visiting scientist * -! * from JMA. * -! * Major modifications: 2002, 2009, 2012 * -! * jim.purser@noaa.gov * -! * * -! ********************************************** -! -! Utility routines for various linear inversions and Cholesky. -! Dependency: modules pkind, pietc -! Originally, these routines were copies of the purely "inversion" members -! of pmat1.f90 (a most extensive collection of matrix routines -- not just -! inversions). As well as having both single and double precision versions -! of each routine, these versions also make provision for a more graceful -! termination in cases where the system matrix is detected to be -! essentially singular (and therefore noninvertible). This provision takes -! the form of an optional "failure flag", FF, which is normally returned -! as .FALSE., but is returned as .TRUE. when inversion fails. -! In Sep 2012, these routines were collected together into pmat.f90 so -! that all the main matrix routines could be in the same library, pmat.a. -! -! DIRECT DEPENDENCIES: -! Modules: pkind, pietc -! -!============================================================================= -module pmat -!============================================================================= -use pkind, only: sp,dp,spc,dpc -use pietc, only: t,f -implicit none -private -public:: ldum,udlmm,inv,L1Lm,LdLm,invl,invu -interface swpvv; module procedure sswpvv,dswpvv,cswpvv; end interface -interface ldum - module procedure sldum,dldum,cldum,sldumf,dldumf,cldumf; end interface -interface udlmm - module procedure sudlmm,dudlmm,cudlmm,sudlmv,dudlmv,cudlmv; end interface -interface inv - module procedure & -sinvmt, dinvmt, cinvmt, slinmmt, dlinmmt, clinmmt, slinmvt, dlinmvt, clinmvt, & -sinvmtf,dinvmtf,cinvmtf,slinmmtf,dlinmmtf,clinmmtf,slinmvtf,dlinmvtf,clinmvtf,& -iinvf - end interface -interface L1Lm; module procedure sL1Lm,dL1Lm,sL1Lmf,dL1Lmf; end interface -interface LdLm; module procedure sLdLm,dLdLm,sLdLmf,dLdLmf; end interface -interface invl; module procedure sinvl,dinvl,slinlv,dlinlv; end interface -interface invu; module procedure sinvu,dinvu,slinuv,dlinuv; end interface - -contains - -!============================================================================= -subroutine sswpvv(d,e)! [swpvv] -!============================================================================= -! Swap vectors -!------------- -real(sp), intent(inout) :: d(:), e(:) -real(sp) :: tv(size(d)) -!============================================================================= -tv = d; d = e; e = tv -end subroutine sswpvv -!============================================================================= -subroutine dswpvv(d,e)! [swpvv] -!============================================================================= -real(dp), intent(inout) :: d(:), e(:) -real(dp) :: tv(size(d)) -!============================================================================= -tv = d; d = e; e = tv -end subroutine dswpvv -!============================================================================= -subroutine cswpvv(d,e)! [swpvv] -!============================================================================= -complex(dpc),intent(inout) :: d(:), e(:) -complex(dpc) :: tv(size(d)) -!============================================================================= -tv = d; d = e; e = tv -end subroutine cswpvv - -!============================================================================= -subroutine sinvmt(a)! [inv] -!============================================================================= -real(sp),dimension(:,:),intent(INOUT):: a -logical :: ff -call sinvmtf(a,ff) -if(ff)stop 'In sinvmt; Unable to invert matrix' -end subroutine sinvmt -!============================================================================= -subroutine dinvmt(a)! [inv] -!============================================================================= -real(dp),dimension(:,:),intent(inout):: a -logical :: ff -call dinvmtf(a,ff) -if(ff)stop 'In dinvmt; Unable to invert matrix' -end subroutine dinvmt -!============================================================================= -subroutine cinvmt(a)! [inv] -!============================================================================= -complex(dpc),dimension(:,:),intent(inout):: a -logical :: ff -call cinvmtf(a,ff) -if(ff)stop 'In cinvmt; Unable to invert matrix' -end subroutine cinvmt -!============================================================================= -subroutine sinvmtf(a,ff)! [inv] -!============================================================================= -! Invert matrix (or flag if can't) -!---------------- -real(sp),dimension(:,:),intent(inout):: a -logical, intent( out):: ff -integer :: m,i,j,jp,l -real(sp) :: d -integer,dimension(size(a,1)) :: ipiv -!============================================================================= -m=size(a,1) -if(m /= size(a,2))stop 'In sinvmtf; matrix passed to sinvmtf is not square' -! Perform a pivoted L-D-U decomposition on matrix a: -call sldumf(a,ipiv,d,ff) -if(ff)then - print '(" In sinvmtf; failed call to sldumf")' - return -endif - -! Invert upper triangular portion U in place: -do i=1,m; a(i,i)=1./a(i,i); enddo -do i=1,m-1 - do j=i+1,m; a(i,j)=-a(j,j)*dot_product(a(i:j-1,j),a(i,i:j-1)); enddo -enddo - -! Invert lower triangular portion L in place: -do j=1,m-1; jp=j+1 - do i=jp,m; a(i,j)=-a(i,j)-dot_product(a(jp:i-1,j),a(i,jp:i-1)); enddo -enddo - -! Form the product of U**-1 and L**-1 in place -do j=1,m-1; jp=j+1 - do i=1,j; a(i,j)=a(i,j)+dot_product(a(jp:m,j),a(i,jp:m)); enddo - do i=jp,m; a(i,j)=dot_product(a(i:m,j),a(i,i:m)); enddo -enddo - -! Permute columns according to ipiv -do j=m-1,1,-1; l=ipiv(j); call sswpvv(a(:,j),a(:,l)); enddo -end subroutine sinvmtf -!============================================================================= -subroutine dinvmtf(a,ff)! [inv] -!============================================================================= -real(DP),dimension(:,:),intent(INOUT):: a -logical, intent( OUT):: ff -integer :: m,i,j,jp,l -real(DP) :: d -integer, dimension(size(a,1)) :: ipiv -!============================================================================= -m=size(a,1) -if(m /= size(a,2))stop 'In inv; matrix passed to dinvmtf is not square' -! Perform a pivoted L-D-U decomposition on matrix a: -call dldumf(a,ipiv,d,ff) -if(ff)then - print '(" In dinvmtf; failed call to dldumf")' - return -endif - -! Invert upper triangular portion U in place: -do i=1,m; a(i,i)=1/a(i,i); enddo -do i=1,m-1 - do j=i+1,m; a(i,j)=-a(j,j)*dot_product(a(i:j-1,j),a(i,i:j-1)); enddo -enddo - -! Invert lower triangular portion L in place: -do j=1,m-1; jp=j+1 - do i=jp,m; a(i,j)=-a(i,j)-dot_product(a(jp:i-1,j),a(i,jp:i-1)); enddo -enddo - -! Form the product of U**-1 and L**-1 in place -do j=1,m-1; jp=j+1 - do i=1,j; a(i,j)=a(i,j)+dot_product(a(jp:m,j),a(i,jp:m)); enddo - do i=jp,m; a(i,j)=dot_product(a(i:m,j),a(i,i:m)); enddo -enddo - -! Permute columns according to ipiv -do j=m-1,1,-1; l=ipiv(j); call dswpvv(a(:,j),a(:,l)); enddo -end subroutine dinvmtf -!============================================================================= -subroutine cinvmtf(a,ff)! [inv] -!============================================================================= -complex(dpc),dimension(:,:),intent(INOUT):: a -logical, intent( OUT):: ff -integer :: m,i,j,jp,l -complex(dpc) :: d -integer, dimension(size(a,1)) :: ipiv -!============================================================================= -m=size(a,1) -if(m /= size(a,2))stop 'In inv; matrix passed to cinvmtf is not square' -! Perform a pivoted L-D-U decomposition on matrix a: -call cldumf(a,ipiv,d,ff) -if(ff)then - print '(" In cinvmtf; failed call to cldumf")' - return -endif - -! Invert upper triangular portion U in place: -do i=1,m; a(i,i)=1/a(i,i); enddo -do i=1,m-1 - do j=i+1,m; a(i,j)=-a(j,j)*sum(a(i:j-1,j)*a(i,i:j-1)); enddo -enddo - -! Invert lower triangular portion L in place: -do j=1,m-1; jp=j+1 - do i=jp,m; a(i,j)=-a(i,j)-sum(a(jp:i-1,j)*a(i,jp:i-1)); enddo -enddo - -! Form the product of U**-1 and L**-1 in place -do j=1,m-1; jp=j+1 - do i=1,j; a(i,j)=a(i,j)+sum(a(jp:m,j)*a(i,jp:m)); enddo - do i=jp,m; a(i,j)=sum(a(i:m,j)*a(i,i:m)); enddo -enddo - -! Permute columns according to ipiv -do j=m-1,1,-1; l=ipiv(j); call cswpvv(a(:,j),a(:,l)); enddo -end subroutine cinvmtf - -!============================================================================= -subroutine slinmmt(a,b)! [inv] -!============================================================================= -real(sp),dimension(:,:),intent(inout):: a,b -logical :: ff -call slinmmtf(a,b,ff) -if(ff)stop 'In slinmmt; unable to invert linear system' -end subroutine slinmmt -!============================================================================= -subroutine dlinmmt(a,b)! [inv] -!============================================================================= -real(dp),dimension(:,:),intent(inout):: a,b -logical :: ff -call dlinmmtf(a,b,ff) -if(ff)stop 'In dlinmmt; unable to invert linear system' -end subroutine dlinmmt -!============================================================================= -subroutine clinmmt(a,b)! [inv] -!============================================================================= -complex(dpc),dimension(:,:),intent(inout):: a,b -logical :: ff -call clinmmtf(a,b,ff) -if(ff)stop 'In clinmmt; unable to invert linear system' -end subroutine clinmmt -!============================================================================= -subroutine slinmmtf(a,b,ff)! [inv] -!============================================================================= -real(SP), dimension(:,:),intent(INOUT):: a,b -logical, intent( OUT):: ff -integer,dimension(size(a,1)) :: ipiv -integer :: m -real(sp) :: d -!============================================================================= -m=size(a,1) -if(m /= size(a,2))stop 'In inv; matrix passed to slinmmtf is not square' -if(m /= size(b,1))& - stop 'In inv; matrix and vectors in slinmmtf have unmatched sizes' -call sldumf(a,ipiv,d,ff) -if(ff)then - print '("In slinmmtf; failed call to sldumf")' - return -endif -call sudlmm(a,b,ipiv) -end subroutine slinmmtf -!============================================================================= -subroutine dlinmmtf(a,b,ff)! [inv] -!============================================================================= -real(dp),dimension(:,:), intent(inout):: a,b -logical, intent( out):: ff -integer, dimension(size(a,1)) :: ipiv -integer :: m -real(dp) :: d -!============================================================================= -m=size(a,1) -if(m /= size(a,2))stop 'In inv; matrix passed to dlinmmtf is not square' -if(m /= size(b,1))& - stop 'In inv; matrix and vectors in dlinmmtf have unmatched sizes' -call dldumf(a,ipiv,d,ff) -if(ff)then - print '("In dlinmmtf; failed call to dldumf")' - return -endif -call dudlmm(a,b,ipiv) -end subroutine dlinmmtf -!============================================================================= -subroutine clinmmtf(a,b,ff)! [inv] -!============================================================================= -complex(dpc),dimension(:,:),intent(INOUT):: a,b -logical, intent( OUT):: ff -integer, dimension(size(a,1)) :: ipiv -integer :: m -complex(dpc) :: d -!============================================================================= -m=size(a,1) -if(m /= size(a,2))stop 'In inv; matrix passed to dlinmmtf is not square' -if(m /= size(b,1))& - stop 'In inv; matrix and vectors in dlinmmtf have unmatched sizes' -call cldumf(a,ipiv,d,ff) -if(ff)then - print '("In clinmmtf; failed call to cldumf")' - return -endif -call cudlmm(a,b,ipiv) -end subroutine clinmmtf - -!============================================================================= -subroutine slinmvt(a,b)! [inv] -!============================================================================= -real(sp), dimension(:,:),intent(inout):: a -real(sp), dimension(:), intent(inout):: b -logical :: ff -call slinmvtf(a,b,ff) -if(ff)stop 'In slinmvt; matrix singular, unable to continue' -end subroutine slinmvt -!============================================================================= -subroutine dlinmvt(a,b)! [inv] -!============================================================================= -real(dp), dimension(:,:),intent(inout):: a -real(dp), dimension(:), intent(inout):: b -logical :: ff -call dlinmvtf(a,b,ff) -if(ff)stop 'In dlinmvt; matrix singular, unable to continue' -end subroutine dlinmvt -!============================================================================= -subroutine clinmvt(a,b)! [inv] -!============================================================================= -complex(dpc), dimension(:,:),intent(inout):: a -complex(dpc), dimension(:), intent(inout):: b -logical :: ff -call clinmvtf(a,b,ff) -if(ff)stop 'In clinmvt; matrix singular, unable to continue' -end subroutine clinmvt -!============================================================================= -subroutine slinmvtf(a,b,ff)! [inv] -!============================================================================= -real(sp),dimension(:,:),intent(inout):: a -real(sp),dimension(:), intent(inout):: b -logical, intent( out):: ff -integer,dimension(size(a,1)) :: ipiv -real(sp) :: d -!============================================================================= -if(size(a,1) /= size(a,2).or. size(a,1) /= size(b))& - stop 'In inv; In slinmvtf; incompatible array dimensions' -call sldumf(a,ipiv,d,ff) -if(ff)then - print '("In slinmvtf; failed call to sldumf")' - return -endif -call sudlmv(a,b,ipiv) -end subroutine slinmvtf -!============================================================================= -subroutine dlinmvtf(a,b,ff)! [inv] -!============================================================================= -real(dp),dimension(:,:),intent(inout):: a -real(dp),dimension(:), intent(inout):: b -logical, intent( out):: ff -integer, dimension(size(a,1)) :: ipiv -real(dp) :: d -!============================================================================= -if(size(a,1) /= size(a,2).or. size(a,1) /= size(b))& - stop 'In inv; incompatible array dimensions passed to dlinmvtf' -call dldumf(a,ipiv,d,ff) -if(ff)then - print '("In dlinmvtf; failed call to dldumf")' - return -endif -call dudlmv(a,b,ipiv) -end subroutine dlinmvtf -!============================================================================= -subroutine clinmvtf(a,b,ff)! [inv] -!============================================================================= -complex(dpc),dimension(:,:),intent(inout):: a -complex(dpc),dimension(:), intent(inout):: b -logical, intent( out):: ff -integer, dimension(size(a,1)) :: ipiv -complex(dpc) :: d -!============================================================================= -if(size(a,1) /= size(a,2).or. size(a,1) /= size(b))& - stop 'In inv; incompatible array dimensions passed to clinmvtf' -call cldumf(a,ipiv,d,ff) -if(ff)then - print '("In clinmvtf; failed call to cldumf")' - return -endif -call cudlmv(a,b,ipiv) -end subroutine clinmvtf - -!============================================================================= -subroutine iinvf(imat,ff)! [inv] -!============================================================================= -! Invert integer square array, imat, if possible, but flag ff=.true. -! if not possible. (Determinant of imat must be +1 or -1 -!============================================================================= -integer,dimension(:,:),intent(INOUT):: imat -logical, intent( OUT):: ff -!----------------------------------------------------------------------------- -real(dp),parameter :: eps=1.e-10_dp -real(dp),dimension(size(imat,1),size(imat,1)):: dmat -integer :: m,i,j -!============================================================================= -m=size(imat,1) -if(m /= size(imat,2))stop 'In inv; matrix passed to iinvf is not square' -dmat=imat; call inv(dmat,ff) -if(.not.ff)then - do j=1,m - do i=1,m - imat(i,j)=nint(dmat(i,j)); if(abs(dmat(i,j)-imat(i,j))>eps)ff=t - enddo - enddo -endif -end subroutine iinvf - -!============================================================================= -subroutine sldum(a,ipiv,d)! [ldum] -!============================================================================= -real(sp),intent(inout) :: a(:,:) -real(sp),intent(out ) :: d -integer, intent(out ) :: ipiv(:) -logical :: ff -call sldumf(a,ipiv,d,ff) -if(ff)stop 'In sldum; matrix singular, unable to continue' -end subroutine sldum -!============================================================================= -subroutine dldum(a,ipiv,d)! [ldum] -!============================================================================= -real(dp),intent(inout) :: a(:,:) -real(dp),intent(out ) :: d -integer, intent(out ) :: ipiv(:) -logical:: ff -call dldumf(a,ipiv,d,ff) -if(ff)stop 'In dldum; matrix singular, unable to continue' -end subroutine dldum -!============================================================================= -subroutine cldum(a,ipiv,d)! [ldum] -!============================================================================= -complex(dpc),intent(inout) :: a(:,:) -complex(dpc),intent(out ) :: d -integer, intent(out ) :: ipiv(:) -logical:: ff -call cldumf(a,ipiv,d,ff) -if(ff)stop 'In cldum; matrix singular, unable to continue' -end subroutine cldum -!============================================================================= -subroutine sldumf(a,ipiv,d,ff)! [ldum] -!============================================================================= -! R.J.Purser, NCEP, Washington D.C. 1996 -! SUBROUTINE LDUM -! perform l-d-u decomposition of square matrix a in place with -! pivoting. -! -! <-> a square matrix to be factorized -! <-- ipiv array encoding the pivoting sequence -! <-- d indicator for possible sign change of determinant -! <-- ff: failure flag, set to .true. when determinant of a vanishes. -!============================================================================= -real(SP),intent(INOUT) :: a(:,:) -real(SP),intent(OUT ) :: d -integer, intent(OUT ) :: ipiv(:) -logical, intent(OUT ) :: ff -integer :: m,i, j, jp, ibig, jm -real(SP) :: s(size(a,1)), aam, aa, abig, ajj, ajji, aij -!============================================================================= -ff=f -m=size(a,1) -do i=1,m - aam=0 - do j=1,m - aa=abs(a(i,j)) - if(aa > aam)aam=aa - enddo - if(aam == 0)then - print '("In sldumf; row ",i6," of matrix vanishes")',i - ff=t - return - endif - s(i)=1/aam -enddo -d=1. -ipiv(m)=m -do j=1,m-1 - jp=j+1 - abig=s(j)*abs(a(j,j)) - ibig=j - do i=jp,m - aa=s(i)*abs(a(i,j)) - if(aa > abig)then - ibig=i - abig=aa - endif - enddo -! swap rows, recording changed sign of determinant - ipiv(j)=ibig - if(ibig /= j)then - d=-d - call sswpvv(a(j,:),a(ibig,:)) - s(ibig)=s(j) - endif - ajj=a(j,j) - if(ajj == 0)then - jm=j-1 - print '(" failure in sldumf:"/" matrix singular, rank=",i3)',jm - ff=t - return - endif - ajji=1/ajj - do i=jp,m - aij=ajji*a(i,j) - a(i,j)=aij - a(i,jp:m) = a(i,jp:m) - aij*a(j,jp:m) - enddo -enddo -end subroutine sldumf -!============================================================================= -subroutine DLDUMf(A,IPIV,D,ff)! [ldum] -!============================================================================= -real(DP), intent(INOUT) :: a(:,:) -real(DP), intent(OUT ) :: d -integer, intent(OUT ) :: ipiv(:) -logical, intent(OUT ) :: ff -integer :: m,i, j, jp, ibig, jm -real(DP) :: s(size(a,1)), aam, aa, abig, ajj, ajji, aij -!============================================================================= -ff=f -m=size(a,1) -do i=1,m - aam=0 - do j=1,m - aa=abs(a(i,j)) - if(aa > aam)aam=aa - enddo - if(aam == 0)then - print '("In dldumf; row ",i6," of matrix vanishes")',i - ff=t - return - endif - s(i)=1/aam -enddo -d=1. -ipiv(m)=m -do j=1,m-1 - jp=j+1 - abig=s(j)*abs(a(j,j)) - ibig=j - do i=jp,m - aa=s(i)*abs(a(i,j)) - if(aa > abig)then - ibig=i - abig=aa - endif - enddo -! swap rows, recording changed sign of determinant - ipiv(j)=ibig - if(ibig /= j)then - d=-d - call dswpvv(a(j,:),a(ibig,:)) - s(ibig)=s(j) - endif - ajj=a(j,j) - if(ajj == 0)then - jm=j-1 - print '(" Failure in dldumf:"/" matrix singular, rank=",i3)',jm - ff=t - return - endif - ajji=1/ajj - do i=jp,m - aij=ajji*a(i,j) - a(i,j)=aij - a(i,jp:m) = a(i,jp:m) - aij*a(j,jp:m) - enddo -enddo -end subroutine DLDUMf -!============================================================================= -subroutine cldumf(a,ipiv,d,ff)! [ldum] -!============================================================================= -use pietc, only: c0 -complex(dpc), intent(INOUT) :: a(:,:) -complex(dpc), intent(OUT ) :: d -integer, intent(OUT ) :: ipiv(:) -logical, intent(OUT ) :: ff -integer :: m,i, j, jp, ibig, jm -complex(dpc) :: ajj, ajji, aij -real(dp) :: aam,aa,abig -real(dp),dimension(size(a,1)):: s -!============================================================================= -ff=f -m=size(a,1) -do i=1,m - aam=0 - do j=1,m - aa=abs(a(i,j)) - if(aa > aam)aam=aa - enddo - if(aam == 0)then - print '("In cldumf; row ",i6," of matrix vanishes")',i - ff=t - return - endif - s(i)=1/aam -enddo -d=1. -ipiv(m)=m -do j=1,m-1 - jp=j+1 - abig=s(j)*abs(a(j,j)) - ibig=j - do i=jp,m - aa=s(i)*abs(a(i,j)) - if(aa > abig)then - ibig=i - abig=aa - endif - enddo -! swap rows, recording changed sign of determinant - ipiv(j)=ibig - if(ibig /= j)then - d=-d - call cswpvv(a(j,:),a(ibig,:)) - s(ibig)=s(j) - endif - ajj=a(j,j) - if(ajj == c0)then - jm=j-1 - print '(" Failure in cldumf:"/" matrix singular, rank=",i3)',jm - ff=t - return - endif - ajji=1/ajj - do i=jp,m - aij=ajji*a(i,j) - a(i,j)=aij - a(i,jp:m) = a(i,jp:m) - aij*a(j,jp:m) - enddo -enddo -end subroutine cldumf - -!============================================================================= -subroutine sudlmm(a,b,ipiv)! [udlmm] -!============================================================================= -! R.J.Purser, National Meteorological Center, Washington D.C. 1993 -! SUBROUTINE UDLMM -! use l-u factors in A to back-substitute for several rhs in B, using ipiv to -! define the pivoting permutation used in the l-u decomposition. -! -! --> A L-D-U factorization of linear system matrux -! <-> B rt-hand-sides vectors on input, corresponding solutions on return -! --> IPIV array encoding the pivoting sequence -!============================================================================= -integer, dimension(:), intent(in) :: ipiv -real(sp),dimension(:,:),intent(in) :: a -real(sp),dimension(:,:),intent(inout) :: b -integer :: m,i, k, l -real(sp) :: s,aiii -!============================================================================= -m=size(a,1) -do k=1,size(b,2) !loop over columns of b - do i=1,m - l=ipiv(i) - s=b(l,k) - b(l,k)=b(i,k) - s = s - sum(b(1:i-1,k)*a(i,1:i-1)) - b(i,k)=s - enddo - b(m,k)=b(m,k)/a(m,m) - do i=m-1,1,-1 - aiii=1/a(i,i) - b(i,k) = b(i,k) - sum(b(i+1:m,k)*a(i,i+1:m)) - b(i,k)=b(i,k)*aiii - enddo -enddo -end subroutine sudlmm -!============================================================================= -subroutine dudlmm(a,b,ipiv)! [udlmm] -!============================================================================= -integer, dimension(:), intent(in ) :: ipiv -real(dp), dimension(:,:),intent(in ) :: a -real(dp), dimension(:,:),intent(inout) :: b -integer :: m,i, k, l -real(dp) :: s,aiii -!============================================================================= -m=size(a,1) -do k=1, size(b,2)!loop over columns of b - do i=1,m - l=ipiv(i) - s=b(l,k) - b(l,k)=b(i,k) - s = s - sum(b(1:i-1,k)*a(i,1:i-1)) - b(i,k)=s - enddo - b(m,k)=b(m,k)/a(m,m) - do i=m-1,1,-1 - aiii=1/a(i,i) - b(i,k) = b(i,k) - sum(b(i+1:m,k)*a(i,i+1:m)) - b(i,k)=b(i,k)*aiii - enddo -enddo -end subroutine dudlmm -!============================================================================= -subroutine cudlmm(a,b,ipiv)! [udlmm] -!============================================================================= -integer, dimension(:), intent(in ) :: ipiv -complex(dpc),dimension(:,:),intent(in ) :: a -complex(dpc),dimension(:,:),intent(inout) :: b -integer :: m,i, k, l -complex(dpc) :: s,aiii -!============================================================================= -m=size(a,1) -do k=1, size(b,2)!loop over columns of b - do i=1,m - l=ipiv(i) - s=b(l,k) - b(l,k)=b(i,k) - s = s - sum(b(1:i-1,k)*a(i,1:i-1)) - b(i,k)=s - enddo - b(m,k)=b(m,k)/a(m,m) - do i=m-1,1,-1 - aiii=1/a(i,i) - b(i,k) = b(i,k) - sum(b(i+1:m,k)*a(i,i+1:m)) - b(i,k)=b(i,k)*aiii - enddo -enddo -end subroutine cudlmm - -!============================================================================= -subroutine sudlmv(a,b,ipiv)! [udlmv] -!============================================================================= -! R.J.Purser, National Meteorological Center, Washington D.C. 1993 -! SUBROUTINE UDLMV -! use l-u factors in A to back-substitute for 1 rhs in B, using ipiv to -! define the pivoting permutation used in the l-u decomposition. -! -! --> A L-D-U factorization of linear system matrix -! <-> B right-hand-side vector on input, corresponding solution on return -! --> IPIV array encoding the pivoting sequence -!============================================================================= -integer, dimension(:), intent(in) :: ipiv -real(sp),dimension(:,:),intent(in) :: a -real(sp),dimension(:), intent(inout) :: b -integer :: m,i, l -real(sp) :: s,aiii -!============================================================================= -m=size(a,1) -do i=1,m - l=ipiv(i) - s=b(l) - b(l)=b(i) - s = s - sum(b(1:i-1)*a(i,1:i-1)) - b(i)=s -enddo -b(m)=b(m)/a(m,m) -do i=m-1,1,-1 - aiii=1/a(i,i) - b(i) = b(i) - sum(b(i+1:m)*a(i,i+1:m)) - b(i)=b(i)*aiii -enddo -end subroutine sudlmv -!============================================================================= -subroutine dudlmv(a,b,ipiv)! [udlmv] -!============================================================================= -integer, dimension(:), intent(in ) :: ipiv(:) -real(dp), dimension(:,:),intent(in ) :: a(:,:) -real(dp), dimension(:), intent(inout) :: b(:) -integer :: m,i, l -real(dp) :: s,aiii -!============================================================================= -m=size(a,1) -do i=1,m - l=ipiv(i) - s=b(l) - b(l)=b(i) - s = s - sum(b(1:i-1)*a(i,1:i-1)) - b(i)=s -enddo -b(m)=b(m)/a(m,m) -do i=m-1,1,-1 - aiii=1/a(i,i) - b(i) = b(i) - sum(b(i+1:m)*a(i,i+1:m)) - b(i)=b(i)*aiii -enddo -end subroutine dudlmv -!============================================================================= -subroutine cudlmv(a,b,ipiv)! [udlmv] -!============================================================================= -integer, dimension(:), intent(in ) :: ipiv(:) -complex(dpc),dimension(:,:),intent(in ) :: a(:,:) -complex(dpc),dimension(:), intent(inout) :: b(:) -integer :: m,i, l -complex(dpc) :: s,aiii -!============================================================================= -m=size(a,1) -do i=1,m - l=ipiv(i) - s=b(l) - b(l)=b(i) - s = s - sum(b(1:i-1)*a(i,1:i-1)) - b(i)=s -enddo -b(m)=b(m)/a(m,m) -do i=m-1,1,-1 - aiii=1/a(i,i) - b(i) = b(i) - sum(b(i+1:m)*a(i,i+1:m)) - b(i)=b(i)*aiii -enddo -end subroutine cudlmv - -!============================================================================= -subroutine sl1lm(a,b) ! [l1lm] -!============================================================================= -! Cholesky, M -> L*U, U(i,j)=L(j,i) -!============================================================================= -real(sp), intent(in ) :: a(:,:) -real(sp), intent(inout) :: b(:,:) -!----------------------------------------------------------------------------- -logical:: ff -call sl1lmf(a,b,ff) -if(ff)stop 'In sl1lm; matrix singular, unable to continue' -end subroutine sl1lm -!============================================================================= -subroutine dl1lm(a,b) ! [l1lm] -!============================================================================= -! Cholesky, M -> L*U, U(i,j)=L(j,i) -!============================================================================= -real(dp), intent(in ) :: a(:,:) -real(dp), intent(inout) :: b(:,:) -!----------------------------------------------------------------------------- -logical:: ff -call dl1lmf(a,b,ff) -if(ff)stop 'In dl1lm; matrix singular, unable to continue' -end subroutine dl1lm - -!============================================================================= -subroutine sl1lmf(a,b,ff)! [L1Lm] -!============================================================================= -! Cholesky, M -> L*U, U(i,j)=L(j,i) -!============================================================================= -real(sp), intent(IN ) :: a(:,:) -real(sp), intent(INOUT) :: b(:,:) -logical :: ff -!----------------------------------------------------------------------------- -integer :: m,j, jm, jp, i -real(sp) :: s, bjji -!============================================================================= -m=size(a,1) -ff=f -do j=1,m - jm=j-1 - jp=j+1 - s = a(j,j) - sum(b(j,1:jm)*b(j,1:jm)) - ff=(S <= 0) - if(ff)then - print '("sL1Lmf detects nonpositive a, rank=",i6)',jm - return - endif - b(j,j)=sqrt(s) - bjji=1/b(j,j) - do i=jp,m - s = a(i,j) - sum(b(i,1:jm)*b(j,1:jm)) - b(i,j)=s*bjji - enddo - b(1:jm,j) = 0 -enddo -end subroutine sl1lmf -!============================================================================= -subroutine dl1lmf(a,b,ff) ! [L1Lm] -!============================================================================= -real(dp), intent(IN ) :: a(:,:) -real(dp), intent(INOUT) :: b(:,:) -logical :: ff -!----------------------------------------------------------------------------- -integer :: m,j, jm, jp, i -real(dp) :: s, bjji -!============================================================================= -m=size(a,1) -ff=f -do j=1,m - jm=j-1 - jp=j+1 - s = a(j,j) - sum(b(j,1:jm)*b(j,1:jm)) - ff=(s <= 0) - if(ff)then - print '("dL1LMF detects nonpositive A, rank=",i6)',jm - return - endif - b(j,j)=sqrt(s) - bjji=1/b(j,j) - do i=jp,m - s = a(i,j) - sum(b(i,1:jm)*b(j,1:jm)) - b(i,j)=s*bjji - enddo - b(1:jm,j) = 0 -enddo -return -end subroutine dl1lmf - -!============================================================================= -subroutine sldlm(a,b,d)! [LdLm] -!============================================================================= -! Modified Cholesky decompose Q --> L*D*U, U(i,j)=L(j,i) -!============================================================================= -real(sp), intent(IN ):: a(:,:) -real(sp), intent(INOUT):: b(:,:) -real(sp), intent( OUT):: d(:) -!----------------------------------------------------------------------------- -logical:: ff -call sldlmf(a,b,d,ff) -if(ff)stop 'In sldlm; matrix singular, unable to continue' -end subroutine sldlm -!============================================================================= -subroutine dldlm(a,b,d)! [LdLm] -!============================================================================= -real(dp), intent(IN ):: a(:,:) -real(dp), intent(INOUT):: b(:,:) -real(dp), intent( OUT):: d(:) -!----------------------------------------------------------------------------- -logical:: ff -call dldlmf(a,b,d,ff) -if(ff)stop 'In dldlm; matrix singular, unable to continue' -end subroutine dldlm - -!============================================================================= -subroutine sldlmf(a,b,d,ff) ! [LDLM] -!============================================================================= -! Modified Cholesky decompose Q --> L*D*U -!============================================================================= -real(sp), intent(IN ):: a(:,:) -real(sp), intent(INOUT):: b(:,:) -real(sp), intent( OUT):: d(:) -logical, intent( OUT):: ff -!----------------------------------------------------------------------------- -integer :: m,j, jm, jp, i -real(sp) :: bjji -!============================================================================= -m=size(a,1) -ff=f -do j=1,m - jm=j-1 - jp=j+1 - d(j)=a(j,j) - sum(b(1:jm,j)*b(j,1:jm)) - b(j,j) = 1 - ff=(d(j) == 0) - if(ff)then - print '("In sldlmf; singularity of matrix detected")' - print '("Rank of matrix: ",i6)',jm - return - endif - bjji=1/d(j) - do i=jp,m - b(j,i)=a(i,j) - dot_product(b(1:jm,j),b(i,1:jm)) - b(i,j)=b(j,i)*bjji - enddo - b(1:jm,j)=0 -enddo -end subroutine sldlmf -!============================================================================= -subroutine dldlmf(a,b,d,ff) ! [LDLM] -!============================================================================= -! Modified Cholesky Q --> L*D*U, U(i,j)=L(j,i) -!============================================================================= -real(dp), intent(IN ) :: a(:,:) -real(dp), intent(INOUT) :: b(:,:) -real(dp), intent( OUT) :: d(:) -logical, intent( OUT) :: ff -!----------------------------------------------------------------------------- -integer :: m,j, jm, jp, i -real(dp) :: bjji -!============================================================================= -m=size(a,1) -ff=f -do j=1,m; jm=j-1; jp=j+1 - d(j)=a(j,j) - sum(b(1:jm,j)*b(j,1:jm)) - b(j,j) = 1 - ff=(d(j) == 0) - if(ff)then - print '("In dldlmf; singularity of matrix detected")' - print '("Rank of matrix: ",i6)',jm - return - endif - bjji=1/d(j) - do i=jp,m - b(j,i)=a(i,j) - dot_product(b(1:jm,j),b(i,1:jm)) - b(i,j)=b(j,i)*bjji - enddo - b(1:jm,j)=0 -enddo -end subroutine dldlmf - -!============================================================================== -subroutine sinvu(a)! [invu] -!============================================================================== -! Invert the upper triangular matrix in place by transposing, calling -! invl, and transposing again. -!============================================================================== -real,dimension(:,:),intent(inout):: a -a=transpose(a); call sinvl(a); a=transpose(a) -end subroutine sinvu -!============================================================================== -subroutine dinvu(a)! [invu] -!============================================================================== -real(dp),dimension(:,:),intent(inout):: a -a=transpose(a); call dinvl(a); a=transpose(a) -end subroutine dinvu -!============================================================================== -subroutine sinvl(a)! [invl] -!============================================================================== -! Invert lower triangular matrix in place -!============================================================================== -real(sp), intent(inout) :: a(:,:) -integer :: m,j, i -m=size(a,1) -do j=m,1,-1 - a(1:j-1,j) = 0.0 - a(j,j)=1./a(j,j) - do i=j+1,m - a(i,j)=-a(i,i)*sum(a(j:i-1,j)*a(i,j:i-1)) - enddo -enddo -end subroutine sinvl -!============================================================================== -subroutine dinvl(a)! [invl] -!============================================================================== -real(dp), intent(inout) :: a(:,:) -integer :: m,j, i -m=size(a,1) -do j=m,1,-1 - a(1:j-1,j) = 0.0 - a(j,j)=1./a(j,j) - do i=j+1,m - a(i,j)=-a(i,i)*sum(a(j:i-1,j)*a(i,j:i-1)) - enddo -enddo -end subroutine dinvl - -!============================================================================== -subroutine slinlv(a,u)! [invl] -!============================================================================== -! Solve linear system involving lower triangular system matrix. -!============================================================================== -real, intent(in ) :: a(:,:) -real, intent(inout) :: u(:) -integer :: i -if(size(a,1) /= size(a,2) .or. size(a,1) /= size(u))& - stop 'In slinlv; incompatible array dimensions' -do i=1,size(u); u(i)=(u(i) - sum(u(:i-1)*a(i,:i-1)))/a(i,i); enddo -end subroutine slinlv -!============================================================================== -subroutine dlinlv(a,u)! [invl] -!============================================================================== -real(dp), intent(in ) :: a(:,:) -real(dp), intent(inout) :: u(:) -integer :: i -if(size(a,1) /= size(a,2) .or. size(a,1) /= size(u))& - stop 'In dlinlv; incompatible array dimensions' -do i=1,size(u); u(i)=(u(i) - sum(u(:i-1)*a(i,:i-1)))/a(i,i); enddo -end subroutine dlinlv - -!============================================================================== -subroutine slinuv(a,u)! [invu] -!============================================================================== -! Solve linear system involving upper triangular system matrix. -!============================================================================== -real, intent(in ) :: a(:,:) -real, intent(inout) :: u(:) -integer :: i -if(size(a,1) /= size(a,2) .or. size(a,1) /= size(u))& - stop 'In linuv; incompatible array dimensions' -do i=size(u),1,-1; u(i)=(u(i) - sum(a(i+1:,i)*u(i+1:)))/a(i,i); enddo -end subroutine slinuv -!============================================================================== -subroutine dlinuv(a,u)! [invu] -!============================================================================== -real(dp), intent(in ) :: a(:,:) -real(dp), intent(inout) :: u(:) -integer :: i -if(size(a,1) /= size(a,2) .or. size(a,1) /= size(u))& - stop 'In dlinuv; incompatible array dimensions' -do i=size(u),1,-1; u(i)=(u(i) - sum(a(i+1:,i)*u(i+1:)))/a(i,i); enddo -end subroutine dlinuv - -end module pmat - diff --git a/sorc/regional_grid.fd/pmat4.f90 b/sorc/regional_grid.fd/pmat4.f90 deleted file mode 100644 index 8cb2fcb70..000000000 --- a/sorc/regional_grid.fd/pmat4.f90 +++ /dev/null @@ -1,1924 +0,0 @@ -! -! ********************************************** -! * MODULE pmat4 * -! * R. J. Purser, NOAA/NCEP/EMC Oct 2005 * -! * 18th May 2012 * -! * jim.purser@noaa.gov * -! * * -! ********************************************** -! -! Euclidean geometry, geometric (stereographic) projections, -! related transformations (Mobius). -! Package for handy vector and matrix operations in Euclidean geometry. -! This package is primarily intended for 3D operations and three of the -! functions (Cross_product, Triple_product and Axial) do not possess simple -! generalizations to a generic number N of dimensions. The others, while -! admitting such N-dimensional generalizations, have not all been provided -! with such generic forms here at the time of writing, though some of these -! may be added at a future date. -! -! May 2017: Added routines to facilitate manipulation of 3D rotations, -! their representations by axial vectors, and routines to compute the -! exponentials of matrices (without resort to eigen methods). Also added -! Quaternion and spinor representations of 3D rotations, and their -! conversion routines. -! -! FUNCTION: -! absv: Absolute magnitude of vector as its euclidean length -! Normalized: Normalized version of given real vector -! Orthogonalized: Orthogonalized version of second vector rel. to first unit v. -! Cross_product: Vector cross-product of the given 2 vectors -! Outer_product: outer-product matrix of the given 2 vectors -! Triple_product: Scalar triple product of given 3 vectors -! Det: Determinant of given matrix -! Axial: Convert axial-vector <--> 2-form (antisymmetric matrix) -! Diag: Diagnl of given matrix, or diagonal matrix of given elements -! Trace: Trace of given matrix -! Identity: Identity 3*3 matrix, or identity n*n matrix for a given n -! Sarea: Spherical area subtended by three vectors, or by lat-lon -! increments forming a triangle or quadrilateral -! Huarea: Spherical area subtended by right-angled spherical triangle -! SUBROUTINE: -! Gram: Right-handed orthogonal basis and rank, nrank. The first -! nrank basis vectors span the column range of matrix given, -! OR ("plain" version) simple unpivoted Gram-Schmidt of a -! square matrix. -! -! In addition, we include routines that relate to stereographic projections -! and some associated mobius transformation utilities, since these complex -! operations have a strong geometrical flavor. -! -! DIRECT DEPENDENCIES -! Libraries[their Modules]: pmat[pmat] -! Additional Modules : pkind, pietc -! -!============================================================================ -module pmat4 -!============================================================================ -use pkind, only: sp,dp,dpc -implicit none -private -public:: absv,normalized,orthogonalized, & - cross_product,outer_product,triple_product,det,axial, & - diag,trace,identity,sarea,huarea,dlltoxy, & - normalize,gram,rowops,corral, & - axtoq,qtoax, & - rottoax,axtorot,spintoq,qtospin,rottoq,qtorot,mulqq, & - expmat,zntay,znfun, & - ctoz,ztoc,setmobius, & - mobius,mobiusi - -interface absv; module procedure absv_s,absv_d; end interface -interface normalized;module procedure normalized_s,normalized_d;end interface -interface orthogonalized - module procedure orthogonalized_s,orthogonalized_d; end interface -interface cross_product - module procedure cross_product_s,cross_product_d; end interface -interface outer_product - module procedure outer_product_s,outer_product_d,outer_product_i - end interface -interface triple_product - module procedure triple_product_s,triple_product_d; end interface -interface det; module procedure det_s,det_d,det_i,det_id; end interface -interface axial - module procedure axial3_s,axial3_d,axial33_s,axial33_d; end interface -interface diag - module procedure diagn_s,diagn_d,diagn_i,diagnn_s,diagnn_d,diagnn_i - end interface -interface trace; module procedure trace_s,trace_d,trace_i; end interface -interface identity; module procedure identity_i,identity3_i; end interface -interface huarea; module procedure huarea_s,huarea_d; end interface -interface sarea - module procedure sarea_s,sarea_d,dtarea_s,dtarea_d,dqarea_s,dqarea_d - end interface -interface dlltoxy; module procedure dlltoxy_s,dlltoxy_d; end interface -interface hav; module procedure hav_s, hav_d; end interface -interface normalize;module procedure normalize_s,normalize_d; end interface -interface gram - module procedure gram_s,gram_d,graml_d,plaingram_s,plaingram_d,rowgram - end interface -interface rowops; module procedure rowops; end interface -interface corral; module procedure corral; end interface -interface rottoax; module procedure rottoax; end interface -interface axtorot; module procedure axtorot; end interface -interface spintoq; module procedure spintoq; end interface -interface qtospin; module procedure qtospin; end interface -interface rottoq; module procedure rottoq; end interface -interface qtorot; module procedure qtorot; end interface -interface axtoq; module procedure axtoq; end interface -interface qtoax; module procedure qtoax; end interface -interface setem; module procedure setem; end interface -interface mulqq; module procedure mulqq; end interface -interface expmat; module procedure expmat,expmatd,expmatdd; end interface -interface zntay; module procedure zntay; end interface -interface znfun; module procedure znfun; end interface -interface ctoz; module procedure ctoz; end interface -interface ztoc; module procedure ztoc,ztocd; end interface -interface setmobius;module procedure setmobius,zsetmobius; end interface -interface mobius; module procedure zmobius,cmobius; end interface -interface mobiusi; module procedure zmobiusi; end interface - -contains - -!============================================================================= -function absv_s(a)result(s)! [absv] -!============================================================================= -real(sp),dimension(:),intent(in):: a -real(sp) :: s -s=sqrt(dot_product(a,a)) -end function absv_s -!============================================================================= -function absv_d(a)result(s)! [absv] -!============================================================================= -real(dp),dimension(:),intent(in):: a -real(dp) :: s -s=sqrt(dot_product(a,a)) -end function absv_d - -!============================================================================= -function normalized_s(a)result(b)! [normalized] -!============================================================================= -real(sp),dimension(:),intent(IN):: a -real(sp),dimension(size(a)) :: b -real(sp) :: s -s=absv_s(a); if(s==0)then; b=0;else;b=a/s;endif -end function normalized_s -!============================================================================= -function normalized_d(a)result(b)! [normalized] -!============================================================================= -real(dp),dimension(:),intent(IN):: a -real(dp),dimension(size(a)) :: b -real(dp) :: s -s=absv_d(a); if(s==0)then; b=0;else;b=a/s;endif -end function normalized_d - -!============================================================================= -function orthogonalized_s(u,a)result(b)! [orthogonalized] -!============================================================================= -real(sp),dimension(:),intent(in):: u,a -real(sp),dimension(size(u)) :: b -real(sp) :: s -! Note: this routine assumes u is already normalized -s=dot_product(u,a); b=a-u*s -end function orthogonalized_s -!============================================================================= -function orthogonalized_d(u,a)result(b)! [orthogonalized] -!============================================================================= -real(dp),dimension(:),intent(in):: u,a -real(dp),dimension(size(u)) :: b -real(dp) :: s -! Note: this routine assumes u is already normalized -s=dot_product(u,a); b=a-u*s -end function orthogonalized_d - -!============================================================================= -function cross_product_s(a,b)result(c)! [cross_product] -!============================================================================= -real(sp),dimension(3),intent(in):: a,b -real(sp),dimension(3) :: c -c(1)=a(2)*b(3)-a(3)*b(2); c(2)=a(3)*b(1)-a(1)*b(3); c(3)=a(1)*b(2)-a(2)*b(1) -end function cross_product_s -!============================================================================= -function cross_product_d(a,b)result(c)! [cross_product] -!============================================================================= -real(dp),dimension(3),intent(in):: a,b -real(dp),dimension(3) :: c -c(1)=a(2)*b(3)-a(3)*b(2); c(2)=a(3)*b(1)-a(1)*b(3); c(3)=a(1)*b(2)-a(2)*b(1) -end function cross_product_d - -!============================================================================= -function outer_product_s(a,b)result(c)! [outer_product] -!============================================================================= -real(sp),dimension(:), intent(in ):: a -real(sp),dimension(:), intent(in ):: b -real(sp),DIMENSION(size(a),size(b)):: c -integer :: nb,i -nb=size(b) -do i=1,nb; c(:,i)=a*b(i); enddo -end function outer_product_s -!============================================================================= -function outer_product_d(a,b)result(c)! [outer_product] -!============================================================================= -real(dp),dimension(:), intent(in ):: a -real(dp),dimension(:), intent(in ):: b -real(dp),dimension(size(a),size(b)):: c -integer :: nb,i -nb=size(b) -do i=1,nb; c(:,i)=a*b(i); enddo -end function outer_product_d -!============================================================================= -function outer_product_i(a,b)result(c)! [outer_product] -!============================================================================= -integer,dimension(:), intent(in ):: a -integer,dimension(:), intent(in ):: b -integer,dimension(size(a),size(b)):: c -integer :: nb,i -nb=size(b) -do i=1,nb; c(:,i)=a*b(i); enddo -end function outer_product_i - -!============================================================================= -function triple_product_s(a,b,c)result(tripleproduct)! [triple_product] -!============================================================================= -real(sp),dimension(3),intent(IN ):: a,b,c -real(sp) :: tripleproduct -tripleproduct=dot_product( cross_product(a,b),c ) -end function triple_product_s -!============================================================================= -function triple_product_d(a,b,c)result(tripleproduct)! [triple_product] -!============================================================================= -real(dp),dimension(3),intent(IN ):: a,b,c -real(dp) :: tripleproduct -tripleproduct=dot_product( cross_product(a,b),c ) -end function triple_product_d - -!============================================================================= -function det_s(a)result(det)! [det] -!============================================================================= -real(sp),dimension(:,:),intent(IN ) ::a -real(sp) :: det -real(sp),dimension(size(a,1),size(a,1)):: b -integer :: n,nrank -n=size(a,1) -if(n==3)then - det=triple_product(a(:,1),a(:,2),a(:,3)) -else - call gram(a,b,nrank,det) - if(nranknrank)exit - ab(k:m,k:n)=matmul( transpose(a(:,k:m)),b(:,k:n) ) - ii =maxloc( abs( ab(k:m,k:n)) )+k-1 - val=maxval( abs( ab(k:m,k:n)) ) - if(val<=vcrit)then - nrank=k-1 - exit - endif - i=ii(1) - j=ii(2) - tv=b(:,j) - b(:,j)=-b(:,k) - b(:,k)=tv - tv=a(:,i) - a(:,i)=-a(:,k) - a(:,k)=tv - w(k:n)=matmul( transpose(b(:,k:n)),tv ) - b(:,k)=matmul(b(:,k:n),w(k:n) ) - s=dot_product(b(:,k),b(:,k)) - s=sqrt(s) - if(w(k)<0)s=-s - det=det*s - b(:,k)=b(:,k)/s - do l=k,n - do j=l+1,n - s=dot_product(b(:,l),b(:,j)) - b(:,j)=normalized( b(:,j)-b(:,l)*s ) - enddo - enddo -enddo -end subroutine gram_s - -!============================================================================= -subroutine gram_d(as,b,nrank,det)! [gram] -!============================================================================= -real(dp),dimension(:,:),intent(IN ) :: as -real(dp),dimension(:,:),intent(OUT) :: b -integer, intent(OUT) :: nrank -real(dp), intent(OUT) :: det -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -real(dp),parameter :: crit=1.e-9_dp -real(dp),dimension(size(as,1),size(as,2)):: a -real(dp),dimension(size(as,2),size(as,1)):: ab -real(dp),dimension(size(as,1)) :: tv,w -real(dp) :: val,s,vcrit -integer :: i,j,k,l,m,n -integer,dimension(2) :: ii -!============================================================================= -n=size(as,1) -m=size(as,2) -if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions' -a=as -b=identity(n) -det=1 -val=maxval(abs(a)) -if(val==0)then - nrank=0 - return -endif -vcrit=val*crit -nrank=min(n,m) -do k=1,n - if(k>nrank)exit - ab(k:m,k:n)=matmul( transpose(a(:,k:m)),b(:,k:n) ) - ii =maxloc( abs( ab(k:m,k:n)) )+k-1 - val=maxval( abs( ab(k:m,k:n)) ) - if(val<=vcrit)then - nrank=k-1 - exit - endif - i=ii(1) - j=ii(2) - tv=b(:,j) - b(:,j)=-b(:,k) - b(:,k)=tv - tv=a(:,i) - a(:,i)=-a(:,k) - a(:,k)=tv - w(k:n)=matmul( transpose(b(:,k:n)),tv ) - b(:,k)=matmul(b(:,k:n),w(k:n) ) - s=dot_product(b(:,k),b(:,k)) - s=sqrt(s) - if(w(k)<0)s=-s - det=det*s - b(:,k)=b(:,k)/s - do l=k,n - do j=l+1,n - s=dot_product(b(:,l),b(:,j)) - b(:,j)=normalized( b(:,j)-b(:,l)*s ) - enddo - enddo -enddo -end subroutine gram_d - -!============================================================================= -subroutine graml_d(as,b,nrank,detsign,ldet)! [gram] -!============================================================================= -! A version of gram_d where the determinant information is returned in -! logarithmic form (to avoid overflows for large matrices). When the -! matrix is singular, the "sign" of the determinant, detsign, is returned -! as zero (instead of either +1 or -1) and ldet is then just the log of -! the nonzero factors found by the process. -!============================================================================= -real(dp),dimension(:,:),intent(IN ) :: as -real(dp),dimension(:,:),intent(OUT) :: b -integer, intent(OUT) :: nrank -integer, intent(out) :: detsign -real(dp), intent(OUT) :: ldet -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -real(dp),parameter :: crit=1.e-9_dp -real(dp),dimension(size(as,1),size(as,2)):: a -real(dp),dimension(size(as,2),size(as,1)):: ab -real(dp),dimension(size(as,1)) :: tv,w -real(dp) :: val,s,vcrit -integer :: i,j,k,l,m,n -integer,dimension(2) :: ii -!============================================================================= -detsign=1 -n=size(as,1) -m=size(as,2) -if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions' -a=as -b=identity(n) -!det=1 -ldet=0 -val=maxval(abs(a)) -if(val==0)then - nrank=0 - return -endif -vcrit=val*crit -nrank=min(n,m) -do k=1,n - if(k>nrank)exit - ab(k:m,k:n)=matmul( transpose(a(:,k:m)),b(:,k:n) ) - ii =maxloc( abs( ab(k:m,k:n)) )+k-1 - val=maxval( abs( ab(k:m,k:n)) ) - if(val<=vcrit)then - nrank=k-1 - exit - endif - i=ii(1) - j=ii(2) - tv=b(:,j) - b(:,j)=-b(:,k) - b(:,k)=tv - tv=a(:,i) - a(:,i)=-a(:,k) - a(:,k)=tv - w(k:n)=matmul( transpose(b(:,k:n)),tv ) - b(:,k)=matmul(b(:,k:n),w(k:n) ) - s=dot_product(b(:,k),b(:,k)) - s=sqrt(s) - if(w(k)<0)s=-s - if(s<0)then - ldet=ldet+log(-s) - detsign=-detsign - elseif(s>0)then - ldet=ldet+log(s) - else - detsign=0 - endif - -! det=det*s - b(:,k)=b(:,k)/s - do l=k,n - do j=l+1,n - s=dot_product(b(:,l),b(:,j)) - b(:,j)=normalized( b(:,j)-b(:,l)*s ) - enddo - enddo -enddo -end subroutine graml_d - -!============================================================================= -subroutine plaingram_s(b,nrank)! [gram] -!============================================================================= -! A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only. -real(sp),dimension(:,:),intent(INOUT) :: b -integer, intent( OUT) :: nrank -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -real(sp),parameter :: crit=1.e-5_sp -real(sp) :: val,vcrit -integer :: j,k,n -!============================================================================= -n=size(b,1); if(n/=size(b,2))stop 'In gram; matrix needs to be square' -val=maxval(abs(b)) -nrank=0 -if(val==0)then - b=0 - return -endif -vcrit=val*crit -do k=1,n - val=sqrt(dot_product(b(:,k),b(:,k))) - if(val<=vcrit)then - b(:,k:n)=0 - return - endif - b(:,k)=b(:,k)/val - nrank=k - do j=k+1,n - b(:,j)=b(:,j)-b(:,k)*dot_product(b(:,k),b(:,j)) - enddo -enddo -end subroutine plaingram_s - -!============================================================================= -subroutine plaingram_d(b,nrank)! [gram] -!============================================================================= -! A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only. -real(dp),dimension(:,:),intent(INOUT) :: b -integer, intent( OUT) :: nrank -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -real(dp),parameter :: crit=1.e-9_dp -real(dp) :: val,vcrit -integer :: j,k,n -!============================================================================= -n=size(b,1); if(n/=size(b,2))stop 'In gram; matrix needs to be square' -val=maxval(abs(b)) -nrank=0 -if(val==0)then - b=0 - return -endif -vcrit=val*crit -do k=1,n - val=sqrt(dot_product(b(:,k),b(:,k))) - if(val<=vcrit)then - b(:,k:n)=0 - return - endif - b(:,k)=b(:,k)/val - nrank=k - do j=k+1,n - b(:,j)=b(:,j)-b(:,k)*dot_product(b(:,k),b(:,j)) - enddo -enddo -end subroutine plaingram_d - -!============================================================================= -subroutine rowgram(m,n,a,ipiv,tt,b,rank)! [gram] -!============================================================================= -! Without changing (tall) rectangular input matrix a, perform pivoted gram- -! schmidt operations to orthogonalize the rows, until rows that remain become -! negligible. Record the pivoting sequence in ipiv, and the row-normalization -! in tt(j,j) and the row-orthogonalization in tt(i,j), for i>j. Note that -! tt(i,j)=0 for i=n please' -nepss=n*epss -rank=n -aa=a -tt=0 -do ii=1,n - -! At this stage, all rows less than ii are already orthonormalized and are -! orthogonal to all rows at and beyond ii. Find the norms of these lower -! rows and pivot the largest of them into position ii: - maxp=0 - maxi=ii - do i=ii,m - p(i)=dot_product(aa(i,:),aa(i,:)) - if(p(i)>maxp)then - maxp=p(i) - maxi=i - endif - enddo - if(maxpu1(2))then - j=3 -else - j=2 -endif -ss=u1(j) -if(ss==0)stop 'In rotov; invalid rot' -if(j/=2)t1(2,:)=t1(3,:) -t1(2,:)=t1(2,:)/sqrt(ss) - -! Form t1(3,:) as the cross product of t1(1,:) and t1(2,:) -t1(3,1)=t1(1,2)*t1(2,3)-t1(1,3)*t1(2,2) -t1(3,2)=t1(1,3)*t1(2,1)-t1(1,1)*t1(2,3) -t1(3,3)=t1(1,1)*t1(2,2)-t1(1,2)*t1(2,1) - -! Project rot into the frame whose axes are the rows of t1: -t2=matmul(t1,matmul(rot,transpose(t1))) - -! Obtain the rotation angle, gamma, implied by rot, and gammah=gamma/2: -gamma=atan2(t2(2,1),t2(1,1)); gammah=gamma/2 - -! Hence deduce coefficients (in the form of a real 4-vector) of one of the two -! possible equivalent spinors: -s=sin(gammah) -q(0)=cos(gammah) -q(1:3)=t1(3,:)*s -end subroutine rottoq - -!============================================================================== -subroutine qtorot(q,rot)! [qtorot] -!============================================================================== -! Go from quaternion to rotation matrix representations -!============================================================================== -real(dp),dimension(0:3),intent(IN ):: q -real(dp),dimension(3,3),intent(OUT):: rot -!============================================================================= -call setem(q(0),q(1),q(2),q(3),rot) -end subroutine qtorot - -!============================================================================= -subroutine axtoq(v,q)! [axtoq] -!============================================================================= -! Go from an axial 3-vector to its equivalent quaternion -!============================================================================= -real(dp),dimension(3), intent(in ):: v -real(dp),dimension(0:3),intent(out):: q -!----------------------------------------------------------------------------- -real(dp),dimension(3,3):: rot -!============================================================================= -call axtorot(v,rot) -call rottoq(rot,q) -end subroutine axtoq - -!============================================================================= -subroutine qtoax(q,v)! [qtoax] -!============================================================================= -! Go from quaternion to axial 3-vector -!============================================================================= -real(dp),dimension(0:3),intent(in ):: q -real(dp),dimension(3), intent(out):: v -!----------------------------------------------------------------------------- -real(dp),dimension(3,3):: rot -!============================================================================= -call qtorot(q,rot) -call rottoax(rot,v) -end subroutine qtoax - -!============================================================================= -subroutine setem(c,d,e,g,r)! [setem] -!============================================================================= -real(dp), intent(IN ):: c,d,e,g -real(dp),dimension(3,3),intent(OUT):: r -!----------------------------------------------------------------------------- -real(dp):: cc,dd,ee,gg,de,dg,eg,dc,ec,gc -!============================================================================= -cc=c*c; dd=d*d; ee=e*e; gg=g*g -de=d*e; dg=d*g; eg=e*g -dc=d*c; ec=e*c; gc=g*c -r(1,1)=cc+dd-ee-gg; r(2,2)=cc-dd+ee-gg; r(3,3)=cc-dd-ee+gg -r(2,3)=2*(eg-dc); r(3,1)=2*(dg-ec); r(1,2)=2*(de-gc) -r(3,2)=2*(eg+dc); r(1,3)=2*(dg+ec); r(2,1)=2*(de+gc) -end subroutine setem - -!============================================================================= -function mulqq(a,b)result(c)! [mulqq] -!============================================================================= -! Multiply quaternions, a*b, assuming operation performed from right to left -!============================================================================= -real(dp),dimension(0:3),intent(IN ):: a,b -real(dp),dimension(0:3) :: c -!------------------------------------------- -c(0)=a(0)*b(0) -a(1)*b(1) -a(2)*b(2) -a(3)*b(3) -c(1)=a(0)*b(1) +a(1)*b(0) +a(2)*b(3) -a(3)*b(2) -c(2)=a(0)*b(2) +a(2)*b(0) +a(3)*b(1) -a(1)*b(3) -c(3)=a(0)*b(3) +a(3)*b(0) +a(1)*b(2) -a(2)*b(1) -end function mulqq -!============================================================================= -subroutine expmat(n,a,b,detb)! [expmat] -!============================================================================= -! Evaluate the exponential, b, of a matrix, a, of degree n. -! Apply the iterated squaring method, m times, to the approximation to -! exp(a/(2**m)) obtained as a Taylor expansion of degree L -! See Fung, T. C., 2004, Int. J. Numer. Meth. Engng, 59, 1273--1286. -!============================================================================= -use pietc, only: u1,u2,o2 -integer, intent(IN ):: n -real(dp),dimension(n,n),intent(IN ):: a -real(dp),dimension(n,n),intent(OUT):: b -real(dp), intent(OUT):: detb -!----------------------------------------------------------------------------- -integer,parameter :: L=5 -real(dp),dimension(n,n):: c,p -real(dp) :: t -integer :: i,m -!============================================================================= -m=10+log(u1+maxval(abs(a)))/log(u2) -t=o2**m -c=a*t -p=c -b=p -do i=2,L - p=matmul(p,c)/i - b=b+p -enddo -do i=1,m - b=b*2+matmul(b,b) -enddo -do i=1,n - b(i,i)=b(i,i)+1 -enddo -detb=0; do i=1,n; detb=detb+a(i,i); enddo; detb=exp(detb) -end subroutine expmat - -!============================================================================= -subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat] -!============================================================================= -! Like expmat, but for the 1st derivatives also. -!============================================================================= -use pietc, only: u1,u2,o2 -integer, intent(IN ):: n -real(dp),dimension(n,n), intent(IN ):: a -real(dp),dimension(n,n), intent(OUT):: b -real(dp),dimension(n,n,(n*(n+1))/2),intent(OUT):: bd -real(dp), intent(OUT):: detb -real(dp),dimension((n*(n+1))/2), intent(OUT):: detbd -!----------------------------------------------------------------------------- -integer,parameter :: L=5 -real(dp),dimension(n,n) :: c,p -real(dp),dimension(n,n,(n*(n+1))/2):: pd,cd -real(dp) :: t -integer :: i,j,k,m,n1 -!============================================================================= -n1=(n*(n+1))/2 -m=10+log(u1+maxval(abs(a)))/log(u2) -t=o2**m -c=a*t -p=c -pd=0 -do k=1,n - pd(k,k,k)=t -enddo -k=n -do i=1,n-1 - do j=i+1,n - k=k+1 - pd(i,j,k)=t - pd(j,i,k)=t - enddo -enddo -if(k/=n1)stop 'In expmatd; n1 is inconsistent with n' -cd=pd -b=p -bd=pd - -do i=2,L - do k=1,n1 - pd(:,:,k)=(matmul(cd(:,:,k),p)+matmul(c,pd(:,:,k)))/i - enddo - p=matmul(c,p)/i - b=b+p - bd=bd+pd -enddo -do i=1,m - do k=1,n1 - bd(:,:,k)=2*bd(:,:,k)+matmul(bd(:,:,k),b)+matmul(b,bd(:,:,k)) - enddo - b=b*2+matmul(b,b) -enddo -do i=1,n - b(i,i)=b(i,i)+1 -enddo -detb=0; do i=1,n; detb=detb+a(i,i); enddo; detb=exp(detb) -detbd=0; do k=1,n; detbd(k)=detb; enddo -end subroutine expmatd - -!============================================================================= -subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] -!============================================================================= -! Like expmat, but for the 1st and 2nd derivatives also. -!============================================================================= -use pietc, only: u1,u2,o2 -integer, intent(IN ):: n -real(dp),dimension(n,n), intent(IN ):: a -real(dp),dimension(n,n), intent(OUT):: b -real(dp),dimension(n,n,(n*(n+1))/2), intent(OUT):: bd -real(dp),dimension(n,n,(n*(n+1))/2,(n*(n+1))/2),intent(OUT):: bdd -real(dp), intent(OUT):: detb -real(dp),dimension((n*(n+1))/2), intent(OUT):: detbd -real(dp),dimension((n*(n+1))/2,(n*(n+1))/2), intent(OUT):: detbdd -!----------------------------------------------------------------------------- -integer,parameter :: L=5 -real(dp),dimension(n,n) :: c,p -real(dp),dimension(n,n,(n*(n+1))/2) :: pd,cd -real(dp),dimension(n,n,(n*(n+1))/2,(n*(n+1))/2):: pdd,cdd -real(dp) :: t -integer :: i,j,k,ki,kj,m,n1 -!============================================================================= -n1=(n*(n+1))/2 -m=10+log(u1+maxval(abs(a)))/log(u2) -t=o2**m -c=a*t -p=c -pd=0 -pdd=0 -do k=1,n - pd(k,k,k)=t -enddo -k=n -do i=1,n-1 - do j=i+1,n - k=k+1 - pd(i,j,k)=t - pd(j,i,k)=t - enddo -enddo -if(k/=n1)stop 'In expmatd; n1 is inconsistent with n' -cd=pd -cdd=0 -b=p -bd=pd -bdd=0 - -do i=2,L - do ki=1,n1 - do kj=1,n1 - pdd(:,:,ki,kj)=(matmul(cd(:,:,ki),pd(:,:,kj)) & - + matmul(cd(:,:,kj),pd(:,:,ki)) & - + matmul(c,pdd(:,:,ki,kj)))/i - enddo - enddo - do k=1,n1 - pd(:,:,k)=(matmul(cd(:,:,k),p)+matmul(c,pd(:,:,k)))/i - enddo - p=matmul(c,p)/i - b=b+p - bd=bd+pd - bdd=bdd+pdd -enddo -do i=1,m - do ki=1,n1 - do kj=1,n1 - bdd(:,:,ki,kj)=2*bdd(:,:,ki,kj) & - +matmul(bdd(:,:,ki,kj),b) & - +matmul(bd(:,:,ki),bd(:,:,kj)) & - +matmul(bd(:,:,kj),bd(:,:,ki)) & - +matmul(b,bdd(:,:,ki,kj)) - enddo - enddo - do k=1,n1 - bd(:,:,k)=2*bd(:,:,k)+matmul(bd(:,:,k),b)+matmul(b,bd(:,:,k)) - enddo - b=b*2+matmul(b,b) -enddo -do i=1,n - b(i,i)=b(i,i)+1 -enddo -detb=0; do i=1,n; detb=detb+a(i,i); enddo; detb=exp(detb) -detbd=0; do k=1,n; detbd(k)=detb; enddo -detbdd=0; do ki=1,n; do kj=1,n; detbdd(ki,kj)=detb; enddo; enddo -end subroutine expmatdd - -!============================================================================= -subroutine zntay(n,z,zn)! [zntay] -!============================================================================= -integer, intent(IN ):: n -real(dp),intent(IN ):: z -real(dp),intent(OUT):: zn -!----------------------------------------------------------------------------- -integer,parameter :: ni=100 -real(dp),parameter :: eps0=1.e-16 -integer :: i,i2,n2 -real(dp) :: t,eps,z2 -!============================================================================= -z2=z*2 -n2=n*2 -t=1 -do i=1,n - t=t/(i*2-1) -enddo -eps=t*eps0 -zn=t -t=t -do i=1,ni - i2=i*2 - t=t*z2/(i2*(i2+n2-1)) - zn=zn+t - if(abs(t)0)then - zn=cosh(rz2) - znd=sinh(rz2)/rz2 - zndd=(zn-znd)/z2 - znddd=(znd-3*zndd)/z2 - do i=1,n - i2p3=i*2+3 - zn=znd - znd=zndd - zndd=znddd - znddd=(znd-i2p3*zndd)/z2 - enddo - else - zn=cos(rz2) - znd=sin(rz2)/rz2 - zndd=-(zn-znd)/z2 - znddd=-(znd-3*zndd)/z2 - do i=1,n - i2p3=i*2+3 - zn=znd - znd=zndd - zndd=znddd - znddd=-(znd-i2p3*zndd)/z2 - enddo - endif -endif -end subroutine znfun - -!============================================================================= -! Utility code for various Mobius transformations. If aa1,bb1,cc1,dd1 are -! the coefficients for one transformation, and aa2,bb2,cc2,dd2 are the -! coefficients for a second one, then the coefficients for the mapping -! of a test point, zz, by aa1 etc to zw, followed by a mapping of zw, by -! aa2 etc to zv, is equivalent to a single mapping zz-->zv by the transformatn -! with coefficients aa3,bb3,cc3,dd3, such that, as 2*2 complex matrices: -! -! [ aa3, bb3 ] [ aa2, bb2 ] [ aa1, bb1 ] -! [ ] = [ ] * [ ] -! [ cc3, dd3 ] [ cc2, dd2 ] [ cc1, dd1 ] . -! -! Note that the determinant of these matrices is always +1 -! -!============================================================================= -subroutine ctoz(v, z,infz)! [ctoz] -!============================================================================= -real(dp),dimension(3),intent(IN ):: v -complex(dpc), intent(OUT):: z -logical, intent(OUT):: infz -!----------------------------------------------------------------------------- -real(dp),parameter:: zero=0,one=1 -real(dp) :: rr,zzpi -!============================================================================= -infz=.false. -z=cmplx(v(1),v(2),dpc) -if(v(3)>0)then - zzpi=one/(one+v(3)) -else - rr=v(1)**2+v(2)**2 - infz=(rr==zero); if(infz)return ! <- The point is mapped to infinity (90S) - zzpi=(one-v(3))/rr -endif -z=z*zzpi -end subroutine ctoz - -!============================================================================= -subroutine ztoc(z,infz, v)! [ztoc] -!============================================================================= -complex(dpc), intent(IN ):: z -logical, intent(IN ):: infz -real(dp),dimension(3),intent(OUT):: v -!----------------------------------------------------------------------------- -real(dp),parameter:: zero=0,one=1 -real(dp) :: r,q,rs,rsc,rsbi -!============================================================================= -if(infz)then; v=(/zero,zero,-one/); return; endif -r=real(z); q=aimag(z); rs=r*r+q*q -rsc=one-rs -rsbi=one/(one+rs) -v(1)=2*rsbi*r -v(2)=2*rsbi*q -v(3)=rsc*rsbi -end subroutine ztoc - -!============================================================================= -subroutine ztocd(z,infz, v,vd)! [ztoc] -!============================================================================= -! The convention adopted for the complex derivative is that, for a complex -! infinitesimal map displacement, delta_z, the corresponding infinitesimal -! change of cartesian vector position is delta_v given by: -! delta_v = Real(vd*delta_z). -! Thus, by a kind of Cauchy-Riemann relation, Imag(vd)=v CROSS Real(vd). -! THE DERIVATIVE FOR THE IDEAL POINT AT INFINITY HAS NOT BEEN CODED YET!!! -!============================================================================= -complex(dpc), intent(IN ):: z -logical, intent(IN ):: infz -real(dp),dimension(3), intent(OUT):: v -complex(dpc),dimension(3),intent(OUT):: vd -!----------------------------------------------------------------------------- -real(dp),parameter :: zero=0,one=1 -real(dp) :: r,q,rs,rsc,rsbi,rsbis -real(dp),dimension(3):: u1,u2 -integer :: i -!============================================================================= -if(infz)then; v=(/zero,zero,-one/); return; endif -r=real(z); q=aimag(z); rs=r*r+q*q -rsc=one-rs -rsbi=one/(one+rs) -rsbis=rsbi**2 -v(1)=2*rsbi*r -v(2)=2*rsbi*q -v(3)=rsc*rsbi -u1(1)=2*(one+q*q-r*r)*rsbis -u1(2)=-4*r*q*rsbis -u1(3)=-4*r*rsbis -u2=cross_product(v,u1) -do i=1,3 - vd(i)=cmplx(u1(i),-u2(i),dpc) -enddo -end subroutine ztocd - -!============================================================================ -subroutine setmobius(xc0,xc1,xc2, aa,bb,cc,dd)! [setmobius] -!============================================================================ -! Find the Mobius transformation complex coefficients, aa,bb,cc,dd, -! with aa*dd-bb*cc=1, for a standard (north-)polar stereographic transformation -! that takes cartesian point, xc0 to the north pole, xc1 to (lat=0,lon=0), -! xc2 to the south pole (=complex infinity). -!============================================================================ -real(dp),dimension(3),intent(IN ):: xc0,xc1,xc2 -complex(dpc), intent(OUT):: aa,bb,cc,dd -!---------------------------------------------------------------------------- -real(dp),parameter:: zero=0,one=1 -logical :: infz0,infz1,infz2 -complex(dpc) :: z0,z1,z2,z02,z10,z21 -!============================================================================ -call ctoz(xc0,z0,infz0) -call ctoz(xc1,z1,infz1) -call ctoz(xc2,z2,infz2) -z21=z2-z1 -z02=z0-z2 -z10=z1-z0 - -if( (z0==z1.and.infz0.eqv.infz1).or.& - (z1==z2.and.infz1.eqv.infz2).or.& - (z2==z0.and.infz2.eqv.infz0)) & - stop 'In setmobius; anchor points must be distinct' - -if(infz2 .or. (.not.infz0 .and. abs(z0) XE three cartesian components. -! <-- DLAT degrees latitude -! <-- DLON degrees longitude -!============================================================================= -use pietc, only: u0,rtod -real(sp),dimension(3),intent(IN ):: xe -real(sp), intent(OUT):: dlat,dlon -!----------------------------------------------------------------------------- -real(sp) :: r -!============================================================================= -r=sqrt(xe(1)**2+xe(2)**2) -dlat=atan2(xe(3),r)*rtod -if(r==u0)then - dlon=u0 -else - dlon=atan2(xe(2),xe(1))*rtod -endif -end subroutine sctog - -!============================================================================= -subroutine dctog(xe,dlat,dlon)! [ctog] -!============================================================================= -use pietc, only: u0,rtod -real(dp),dimension(3),intent(IN ):: xe -real(dp), intent(OUT):: dlat,dlon -!----------------------------------------------------------------------------- -real(dp) :: r -!============================================================================= -r=sqrt(xe(1)**2+xe(2)**2) -dlat=atan2(xe(3),r)*rtod -if(r==u0)then - dlon=u0 -else - dlon=atan2(xe(2),xe(1))*rtod -endif -end subroutine dctog - -!============================================================================= -subroutine sgtoc(dlat,dlon,xe)! [gtoc] -!============================================================================= -! R.J.Purser, National Meteorological Center, Washington D.C. 1994 -! SUBROUTINE GTOC -! Transform "Geographical" to "Cartesian" coordinates, where the -! geographical coordinates refer to latitude and longitude (degrees) -! and cartesian coordinates are standard earth-centered cartesian -! coordinates: xe(3) pointing north, xe(1) pointing to the 0-meridian. -! --> DLAT degrees latitude -! --> DLON degrees longitude -! <-- XE three cartesian components. -!============================================================================= -use pietc, only: dtor -real(sp), intent(IN ):: dlat,dlon -real(sp),dimension(3),intent(OUT):: xe -!----------------------------------------------------------------------------- -real(sp) :: rlat,rlon,sla,cla,slo,clo -!============================================================================= -rlat=dtor*dlat; rlon=dtor*dlon -sla=sin(rlat); cla=cos(rlat) -slo=sin(rlon); clo=cos(rlon) -xe(1)=cla*clo; xe(2)=cla*slo; xe(3)=sla -end subroutine sgtoc -!============================================================================= -subroutine dgtoc(dlat,dlon,xe)! [gtoc] -!============================================================================= -use pietc, only: dtor -real(dp), intent(IN ):: dlat,dlon -real(dp),dimension(3),intent(OUT):: xe -!----------------------------------------------------------------------------- -real(dp) :: rlat,rlon,sla,cla,slo,clo -!============================================================================= -rlat=dtor*dlat; rlon=dtor*dlon -sla=sin(rlat); cla=cos(rlat) -slo=sin(rlon); clo=cos(rlon) -xe(1)=cla*clo; xe(2)=cla*slo; xe(3)=sla -end subroutine dgtoc -!============================================================================= -subroutine sgtocd(dlat,dlon,xe,dxedlat,dxedlon)! [gtoc] -!============================================================================= -real(sp), intent(IN ):: dlat,dlon -real(sp),dimension(3),intent(OUT):: xe,dxedlat,dxedlon -!----------------------------------------------------------------------------- -real(dp) :: dlat_d,dlon_d -real(dp),dimension(3):: xe_d,dxedlat_d,dxedlon_d -!============================================================================= -dlat_d=dlat; dlon_d=dlon -call dgtocd(dlat_d,dlon_d,xe_d,dxedlat_d,dxedlon_d) -xe =xe_d -dxedlat=dxedlat_d -dxedlon=dxedlon_d -end subroutine sgtocd -!============================================================================= -subroutine dgtocd(dlat,dlon,xe,dxedlat,dxedlon)! [gtoc] -!============================================================================= -use pietc, only: dtor -real(dp), intent(IN ):: dlat,dlon -real(dp),dimension(3),intent(OUT):: xe,dxedlat,dxedlon -!----------------------------------------------------------------------------- -real(dp) :: rlat,rlon,sla,cla,slo,clo -!============================================================================= -rlat=dtor*dlat; rlon=dtor*dlon -sla=sin(rlat); cla=cos(rlat) -slo=sin(rlon); clo=cos(rlon) -xe(1)=cla*clo; xe(2)=cla*slo; xe(3)=sla -dxedlat(1)=-sla*clo; dxedlat(2)=-sla*slo; dxedlat(3)=cla; dxedlat=dxedlat*dtor -dxedlon(1)=-cla*slo; dxedlon(2)= cla*clo; dxedlon(3)=0 ; dxedlon=dxedlon*dtor -end subroutine dgtocd -!============================================================================= -subroutine sgtocdd(dlat,dlon,xe,dxedlat,dxedlon, & - ddxedlatdlat,ddxedlatdlon,ddxedlondlon)! [gtoc] -!============================================================================= -use pietc, only: dtor -real(sp), intent(IN ):: dlat,dlon -real(sp),dimension(3),intent(OUT):: xe,dxedlat,dxedlon, & - ddxedlatdlat,ddxedlatdlon,ddxedlondlon -!----------------------------------------------------------------------------- -real(dp) :: dlat_d,dlon_d -real(dp),dimension(3):: xe_d,dxedlat_d,dxedlon_d, & - ddxedlatdlat_d,ddxedlatdlon_d,ddxedlondlon_d -!============================================================================= -dlat_d=dlat; dlon_d=dlon -call dgtocdd(dlat_d,dlon_d,xe_d,dxedlat_d,dxedlon_d, & - ddxedlatdlat_d,ddxedlatdlon_d,ddxedlondlon_d) -xe =xe_d -dxedlat =dxedlat_d -dxedlon =dxedlon_d -ddxedlatdlat=ddxedlatdlat_d -ddxedlatdlon=ddxedlatdlon_d -ddxedlondlon=ddxedlondlon_d -end subroutine sgtocdd -!============================================================================= -subroutine dgtocdd(dlat,dlon,xe,dxedlat,dxedlon, & - ddxedlatdlat,ddxedlatdlon,ddxedlondlon)! [gtoc] -!============================================================================= -use pietc, only: dtor -real(dp), intent(IN ):: dlat,dlon -real(dp),dimension(3),intent(OUT):: xe,dxedlat,dxedlon, & - ddxedlatdlat,ddxedlatdlon,ddxedlondlon -!----------------------------------------------------------------------------- -real(dp) :: rlat,rlon,sla,cla,slo,clo -!============================================================================= -rlat=dtor*dlat; rlon=dtor*dlon -sla=sin(rlat); cla=cos(rlat) -slo=sin(rlon); clo=cos(rlon) -xe(1)=cla*clo; xe(2)=cla*slo; xe(3)=sla -dxedlat(1)=-sla*clo; dxedlat(2)=-sla*slo; dxedlat(3)=cla; dxedlat=dxedlat*dtor -dxedlon(1)=-cla*slo; dxedlon(2)= cla*clo; dxedlon(3)=0 ; dxedlon=dxedlon*dtor -ddxedlatdlat(1)=-cla*clo -ddxedlatdlat(2)=-cla*slo -ddxedlatdlat(3)=-sla -ddxedlatdlon(1)= sla*slo -ddxedlatdlon(2)=-sla*clo -ddxedlatdlon(3)= 0 -ddxedlondlon(1)=-cla*clo -ddxedlondlon(2)=-cla*slo -ddxedlondlon(3)= 0 -ddxedlatdlat=ddxedlatdlat*dtor**2 -ddxedlatdlon=ddxedlatdlon*dtor**2 -ddxedlondlon=ddxedlondlon*dtor**2 -end subroutine dgtocdd - -!============================================================================== -subroutine sgtoframem(splat,splon,sorth)! [gtoframe] -!============================================================================== -real(sp), intent(in ):: splat,splon -real(sp),dimension(3,3),intent(out):: sorth -!------------------------------------------------------------------------------ -real(dp):: plat,plon -real(dp),dimension(3,3):: orth -!============================================================================== -plat=splat; plon=splon; call gtoframem(plat,plon,orth); sorth=orth -end subroutine sgtoframem -!============================================================================== -subroutine gtoframem(plat,plon,orth)! [gtoframe] -!============================================================================== -! From the degree lat and lo (plat and plon) return the standard orthogonal -! 3D frame at this location as an orthonormal matrix, orth. -!============================================================================== -real(dp), intent(in ):: plat,plon -real(dp),dimension(3,3),intent(out):: orth -!------------------------------------------------------------------------------ -real(dp),dimension(3):: xp,yp,zp -!============================================================================== -call gtoframev(plat,plon, xp,yp,zp) -orth(:,1)=xp; orth(:,2)=yp; orth(:,3)=zp -end subroutine gtoframem -!============================================================================== -subroutine sgtoframev(splat,splon,sxp,syp,szp)! [gtoframe] -!============================================================================== -real(sp), intent(in ):: splat,splon -real(sp),dimension(3),intent(out):: sxp,syp,szp -!------------------------------------------------------------------------------ -real(dp) :: plat,plon -real(dp),dimension(3):: xp,yp,zp -!============================================================================== -plat=splat; plon=splon -call gtoframev(plat,plon, xp,yp,zp) -sxp=xp; syp=yp; szp=zp -end subroutine sgtoframev -!============================================================================== -subroutine gtoframev(plat,plon, xp,yp,zp)! [gtoframe] -!============================================================================== -! Given a geographical point by its degrees lat and lon, plat and plon, -! return its standard orthogonal cartesian frame, {xp,yp,zp} in earth-centered -! coordinates. -!============================================================================= -use pietc, only: u0,u1 -real(dp), intent(in ):: plat,plon -real(dp),dimension(3),intent(out):: xp,yp,zp -!------------------------------------------------------------------------------ -real(dp),dimension(3):: dzpdlat,dzpdlon -!============================================================================== -if(plat==90)then ! is this the north pole? - xp=(/ u0,u1, u0/) ! Assume the limiting case lat-->90 along the 0-meridian - yp=(/-u1,u0, u0/) ! - zp=(/ u0,u0, u1/) -elseif(plat==-90)then - xp=(/ u0,u1, u0/) ! Assume the limiting case lat-->90 along the 0-meridian - yp=(/ u1,u0, u0/) ! - zp=(/ u0,u0,-u1/) -else - call gtoc(plat,plon,zp,dzpdlat,dzpdlon) - xp=dzpdlon/sqrt(dot_product(dzpdlon,dzpdlon)) - yp=dzpdlat/sqrt(dot_product(dzpdlat,dzpdlat)) -endif -end subroutine gtoframev - -!============================================================================== -subroutine sparaframe(sxp,syp,szp, sxv,syv,szv)! [paraframe] -!============================================================================== -real(sp),dimension(3),intent(in ):: sxp,syp,szp, szv -real(sp),dimension(3),intent(out):: sxv,syv -!----------------------------------------------------------------------------- -real(dp),dimension(3):: xp,yp,zp, xv,yv,zv -!============================================================================= -xp=sxp; yp=syp; zp=szp -call paraframe(xp,yp,zp, xv,yv,zv) -sxv=xv; syv=yv -end subroutine sparaframe -!============================================================================== -subroutine paraframe(xp,yp,zp, xv,yv,zv)! [paraframe] -!============================================================================== -! Take a principal reference orthonormal frame, {xp,yp,zp} and a dependent -! point defined by unit vector, zv, and complete the V-frame cartesian -! components, {xv,yv}, that are the result of parallel-transport of {xp,yp} -! along the geodesic between P and V -!============================================================================== -use pmat4, only: cross_product,normalized -real(dp),dimension(3),intent(in ):: xp,yp,zp, zv -real(dp),dimension(3),intent(out):: xv,yv -!------------------------------------------------------------------------------ -real(dp) :: xpofv,ypofv,theta,ctheta,stheta -real(dp),dimension(3):: xq,yq -!============================================================================== -xpofv=dot_product(xp,zv) -ypofv=dot_product(yp,zv) -theta=atan2(ypofv,xpofv); ctheta=cos(theta); stheta=sin(theta) -xq=zv-zp; xq=xq-zv*dot_product(xq,zv); xq=normalized(xq) -yq=cross_product(zv,xq) -xv=xq*ctheta-yq*stheta -yv=xq*stheta+yq*ctheta -end subroutine paraframe - -!============================================================================== -subroutine sframetwist(sxp,syp,szp, sxv,syv,szv, stwist)! [frametwist] -!============================================================================== -real(sp),dimension(3),intent(in ):: sxp,syp,szp, sxv,syv,szv -real(sp), intent(out):: stwist -!------------------------------------------------------------------------------ -real(dp),dimension(3):: xp,yp,zp, xv,yv,zv -real(dp) :: twist -!============================================================================== -xp=sxp;yp=syp; zp=szp; xv=sxv; yv=syv; zv=szv -call frametwist(xp,yp,zp, xv,yv,zv, twist) -stwist=twist -end subroutine sframetwist -!============================================================================== -subroutine frametwist(xp,yp,zp, xv,yv,zv, twist)! [frametwist] -!============================================================================== -! Given a principal cartesian orthonormal frame, {xp,yp,zp} (i.e., at P with -! Earth-centered cartesians, zp), and another similar frame {xv,yv,zv} at V -! with Earth-centered cartesians zv, find the relative rotation angle, "twist" -! by which the frame at V is rotated in the counterclockwise sense relative -! to the parallel-transportation of P's frame to V. -! Note that, by symmetry, transposing P and V leads to the opposite twist. -!============================================================================== -real(dp),dimension(3),intent(in ):: xp,yp,zp, xv,yv,zv -real(dp), intent(out):: twist -!------------------------------------------------------------------------------ -real(dp),dimension(3):: xxv,yyv -real(dp) :: c,s -!============================================================================== -call paraframe(xp,yp,zp, xxv,yyv,zv) -c=dot_product(xv,xxv); s=dot_product(xv,yyv) -twist=atan2(s,c) -end subroutine frametwist - -!============================================================================= -subroutine sctoc(s,xc1,xc2)! [ctoc_schm] -!============================================================================= -! Evaluate schmidt transformation, xc1 --> xc2, with scaling parameter s -!============================================================================= -real(sp), intent(IN ):: s -real(sp),dimension(3),intent(INOUT):: xc1,xc2 -!----------------------------------------------------------------------------- -real(sp) :: x,y,z,a,b,d,e,ab2,aa,bb,di,aapbb,aambb -!============================================================================= -x=xc1(1); y=xc1(2); z=xc1(3) -a=s+1 -b=s-1 -ab2=a*b*2 -aa=a*a -bb=b*b -aapbb=aa+bb -aambb=aa-bb -d=aapbb-ab2*z -e=aapbb*z-ab2 -di=1/d -xc2(1)=(aambb*x)*di -xc2(2)=(aambb*y)*di -xc2(3)=e*di -end subroutine sctoc - -!============================================================================= -subroutine sctocd(s,xc1,xc2,dxc2)! [ctoc_schm] -!============================================================================= -! Evaluate schmidt transformation, xc1 --> xc2, with scaling parameter s, -! and its jacobian, dxc2. -!============================================================================= -real(sp),intent(IN) :: s -real(sp),dimension(3), intent(INOUT):: xc1,xc2 -real(sp),dimension(3,3),intent( OUT):: dxc2 -!----------------------------------------------------------------------------- -real(sp) :: x,y,z,a,b,d,e, & - ab2,aa,bb,di,ddi,aapbb,aambb -!============================================================================= -x=xc1(1); y=xc1(2); z=xc1(3) -a=s+1 -b=s-1 -ab2=a*b*2 -aa=a*a -bb=b*b -aapbb=aa+bb -aambb=aa-bb -d=aapbb-ab2*z -e=aapbb*z-ab2 -di=1/d -xc2(1)=(aambb*x)*di -xc2(2)=(aambb*y)*di -xc2(3)=e*di -ddi=di*di - -dxc2=0 -dxc2(1,1)=aambb*di -dxc2(1,3)=ab2*aambb*x*ddi -dxc2(2,2)=aambb*di -dxc2(2,3)=ab2*aambb*y*ddi -dxc2(3,3)=aapbb*di +ab2*e*ddi -end subroutine sctocd - -!============================================================================= -subroutine sctocdd(s,xc1,xc2,dxc2,ddxc2)! [ctoc_schm] -!============================================================================= -! Evaluate schmidt transformation, xc1 --> xc2, with scaling parameter s, -! its jacobian, dxc2, and its 2nd derivative, ddxc2. -!============================================================================= -real(sp), intent(IN ):: s -real(sp),dimension(3), intent(INOUT):: xc1,xc2 -real(sp),dimension(3,3), intent( OUT):: dxc2 -real(sp),dimension(3,3,3),intent( OUT):: ddxc2 -!----------------------------------------------------------------------------- -real(sp) :: x,y,z,a,b,d,e, & - ab2,aa,bb,di,ddi,dddi, & - aapbb,aambb -!============================================================================= -x=xc1(1); y=xc1(2); z=xc1(3) -a=s+1 -b=s-1 -ab2=a*b*2 -aa=a*a -bb=b*b -aapbb=aa+bb -aambb=aa-bb -d=aapbb-ab2*z -e=aapbb*z-ab2 -di=1/d -xc2(1)=(aambb*x)*di -xc2(2)=(aambb*y)*di -xc2(3)=e*di -ddi=di*di -dddi=ddi*di - -dxc2=0 -dxc2(1,1)=aambb*di -dxc2(1,3)=ab2*aambb*x*ddi -dxc2(2,2)=aambb*di -dxc2(2,3)=ab2*aambb*y*ddi -dxc2(3,3)=aapbb*di +ab2*e*ddi - -ddxc2=0 -ddxc2(1,1,3)=ab2*aambb*ddi -ddxc2(1,3,1)=ddxc2(1,1,3) -ddxc2(1,3,3)=2*ab2**2*aambb*x*dddi -ddxc2(2,2,3)=ab2*aambb*ddi -ddxc2(2,3,2)=ddxc2(2,2,3) -ddxc2(2,3,3)=2*ab2**2*aambb*y*dddi -ddxc2(3,3,3)=2*ab2*(aapbb*ddi+ab2*e*dddi) -end subroutine sctocdd - -!============================================================================= -subroutine dctoc(s,xc1,xc2)! [ctoc_schm] -!============================================================================= -! Evaluate schmidt transformation, xc1 --> xc2, with scaling parameter s -!============================================================================= -real(dp), intent(IN ):: s -real(dp),dimension(3),intent(INOUT):: xc1,xc2 -!----------------------------------------------------------------------------- -real(dp) :: x,y,z,a,b,d,e, & - ab2,aa,bb,di,aapbb,aambb -!============================================================================= -x=xc1(1); y=xc1(2); z=xc1(3) -a=s+1 -b=s-1 -ab2=a*b*2 -aa=a*a -bb=b*b -aapbb=aa+bb -aambb=aa-bb -d=aapbb-ab2*z -e=aapbb*z-ab2 -di=1/d -xc2(1)=(aambb*x)*di -xc2(2)=(aambb*y)*di -xc2(3)=e*di -end subroutine dctoc - -!============================================================================= -subroutine dctocd(s,xc1,xc2,dxc2)! [ctoc_schm] -!============================================================================= -! Evaluate schmidt transformation, xc1 --> xc2, with scaling parameter s, -! and its jacobian, dxc2. -!============================================================================= -real(dp), intent(IN ):: s -real(dp),dimension(3), intent(INOUT):: xc1,xc2 -real(dp),dimension(3,3),intent( OUT):: dxc2 -!----------------------------------------------------------------------------- -real(dp) :: x,y,z,a,b,d,e, & - ab2,aa,bb,di,ddi,aapbb,aambb -!============================================================================= -x=xc1(1); y=xc1(2); z=xc1(3) -a=s+1 -b=s-1 -ab2=a*b*2 -aa=a*a -bb=b*b -aapbb=aa+bb -aambb=aa-bb -d=aapbb-ab2*z -e=aapbb*z-ab2 -di=1/d -xc2(1)=(aambb*x)*di -xc2(2)=(aambb*y)*di -xc2(3)=e*di -ddi=di*di - -dxc2=0 -dxc2(1,1)=aambb*di -dxc2(1,3)=ab2*aambb*x*ddi -dxc2(2,2)=aambb*di -dxc2(2,3)=ab2*aambb*y*ddi -dxc2(3,3)=aapbb*di +ab2*e*ddi -end subroutine dctocd - -!============================================================================= -subroutine dctocdd(s,xc1,xc2,dxc2,ddxc2)! [ctoc_schm] -!============================================================================= -! Evaluate schmidt transformation, xc1 --> xc2, with scaling parameter s, -! its jacobian, dxc2, and its 2nd derivative, ddxc2. -!============================================================================= -real(dp),intent(IN) :: s -real(dp),dimension(3), intent(INOUT):: xc1,xc2 -real(dp),dimension(3,3), intent(OUT ):: dxc2 -real(dp),dimension(3,3,3),intent(OUT ):: ddxc2 -!----------------------------------------------------------------------------- -real(dp) :: x,y,z,a,b,d,e, & - ab2,aa,bb,di,ddi,dddi, & - aapbb,aambb -!============================================================================= -x=xc1(1); y=xc1(2); z=xc1(3) -a=s+1 -b=s-1 -ab2=a*b*2 -aa=a*a -bb=b*b -aapbb=aa+bb -aambb=aa-bb -d=aapbb-ab2*z -e=aapbb*z-ab2 -di=1/d -xc2(1)=(aambb*x)*di -xc2(2)=(aambb*y)*di -xc2(3)=e*di -ddi=di*di -dddi=ddi*di - -dxc2=0 -dxc2(1,1)=aambb*di -dxc2(1,3)=ab2*aambb*x*ddi -dxc2(2,2)=aambb*di -dxc2(2,3)=ab2*aambb*y*ddi -dxc2(3,3)=aapbb*di +ab2*e*ddi - -ddxc2=0 -ddxc2(1,1,3)=ab2*aambb*ddi -ddxc2(1,3,1)=ddxc2(1,1,3) -ddxc2(1,3,3)=2*ab2**2*aambb*x*dddi -ddxc2(2,2,3)=ab2*aambb*ddi -ddxc2(2,3,2)=ddxc2(2,2,3) -ddxc2(2,3,3)=2*ab2**2*aambb*y*dddi -ddxc2(3,3,3)=2*ab2*(aapbb*ddi+ab2*e*dddi) -end subroutine dctocdd - -!============================================================================= -subroutine plrot(rot3,n,x,y,z)! [plrot] -!============================================================================= -! Apply a constant rotation to a three dimensional polyline -!============================================================================= -integer, intent(IN ):: n -real(sp),dimension(3,3),intent(IN ):: rot3 -real(sp),dimension(n), intent(INOUT):: x,y,z -!----------------------------------------------------------------------------- -real(sp),dimension(3) :: t -integer :: k -!============================================================================= -do k=1,n - t(1)=x(k); t(2)=y(k); t(3)=z(k) - t=matmul(rot3,t) - x(k)=t(1); y(k)=t(2); z(k)=t(3) -enddo -end subroutine plrot - -!============================================================================= -subroutine plroti(rot3,n,x,y,z)! [plroti] -!============================================================================= -! Invert the rotation of a three-dimensional polyline -!============================================================================= -integer, intent(IN ):: n -real(sp),dimension(3,3),intent(IN ):: rot3 -real(sp),dimension(n), intent(INOUT):: x,y,z -!----------------------------------------------------------------------------- -real(sp),dimension(3) :: t -integer :: k -!============================================================================= -do k=1,n - t(1)=x(k); t(2)=y(k); t(3)=z(k) - t=matmul(t,rot3) - x(k)=t(1); y(k)=t(2); z(k)=t(3) -enddo -end subroutine plroti - -!============================================================================= -subroutine dplrot(rot3,n,x,y,z)! [plrot] -!============================================================================= -! Apply a constant rotation to a three dimensional polyline -!============================================================================= -integer, intent(IN ):: n -real(dP),dimension(3,3),intent(IN ):: rot3 -real(dP),dimension(n), intent(INOUT):: x,y,z -!----------------------------------------------------------------------------- -real(dP),dimension(3) :: t -integer :: k -!============================================================================= -do k=1,n - t(1)=x(k); t(2)=y(k); t(3)=z(k) - t=matmul(rot3,t) - x(k)=t(1); y(k)=t(2); z(k)=t(3) -enddo -end subroutine dplrot - -!============================================================================= -subroutine dplroti(rot3,n,x,y,z)! [plroti] -!============================================================================= -! Invert the rotation of a three-dimensional polyline -!============================================================================= -integer, intent(IN ):: n -real(dP),dimension(3,3),intent(IN ):: rot3 -real(dP),dimension(n), intent(INOUT):: x,y,z -!----------------------------------------------------------------------------- -real(dP),dimension(3) :: t -integer :: k -!============================================================================= -do k=1,n - t(1)=x(k); t(2)=y(k); t(3)=z(k) - t=matmul(t,rot3) - x(k)=t(1); y(k)=t(2); z(k)=t(3) -enddo -end subroutine dplroti - -!============================================================================= -subroutine plctoc(s,n,x,y,z)! [plctoc] -!============================================================================= -! Perform schmidt transformation with scaling parameter s to a polyline -!============================================================================= -integer, intent(IN ):: n -real(sp), intent(IN ):: s -real(sp),dimension(n),intent(INOUT):: x,y,z -!----------------------------------------------------------------------------- -real(sp) :: a,b,d,e,ab2,aa,bb,di,aapbb,aambb -integer :: i -!============================================================================= -a=s+1 -b=s-1 -ab2=a*b*2 -aa=a*a -bb=b*b -aapbb=aa+bb -aambb=aa-bb -do i=1,n - d=aapbb-ab2*z(i) - e=aapbb*z(i)-ab2 - di=1/d - x(i)=(aambb*x(i))*di - y(i)=(aambb*y(i))*di - z(i)=e*di -enddo -end subroutine plctoc - -end module pmat5 - - - diff --git a/sorc/regional_grid.fd/psym2.f90 b/sorc/regional_grid.fd/psym2.f90 deleted file mode 100644 index 3b6459047..000000000 --- a/sorc/regional_grid.fd/psym2.f90 +++ /dev/null @@ -1,498 +0,0 @@ -! *********************************** -! * module psym2 * -! * R. J. Purser * -! * NOAA/NCEP/EMC September 2018 * -! * jim.purser@noaa.gov * -! *********************************** -! -! A suite of routines to perform the eigen-decomposition of symmetric 2*2 -! matrices and to deliver basic analytic functions, and the derivatives -! of these functions, of such matrices. -! -! DIRECT DEPENDENCIES -! Module: pkind, pietc -! -!============================================================================= -module psym2 -!============================================================================= -use pkind, only: dp -use pietc, only: u0,u1,o2 -implicit none -private -public:: eigensym2,invsym2,sqrtsym2,expsym2,logsym2,id2222 - -real(dp),dimension(2,2,2,2):: id -data id/u1,u0,u0,u0, u0,o2,o2,u0, u0,o2,o2,u0, u0,u0,u0,u1/! Effective identity - -interface eigensym2; module procedure eigensym2,eigensym2d; end interface -interface invsym2; module procedure invsym2,invsym2d; end interface -interface sqrtsym2; module procedure sqrtsym2,sqrtsym2d; end interface -interface sqrtsym2d_e; module procedure sqrtsym2d_e; end interface -interface sqrtsym2d_t; module procedure sqrtsym2d_t; end interface -interface expsym2; module procedure expsym2,expsym2d; end interface -interface expsym2d_e; module procedure expsym2d_e; end interface -interface expsym2d_t; module procedure expsym2d_t; end interface -interface logsym2; module procedure logsym2,logsym2d ; end interface -interface logsym2d_e; module procedure logsym2d_e; end interface -interface logsym2d_t; module procedure logsym2d_t; end interface -interface id2222; module procedure id2222; end interface - -contains - -!============================================================================= -subroutine eigensym2(em,vv,oo)! [eigensym2] -!============================================================================= -! Get the orthogonal eigenvectors, vv, and diagonal matrix of eigenvalues, oo, -! of the symmetric 2*2 matrix, em. -!============================================================================= -real(dp),dimension(2,2),intent(in ):: em -real(dp),dimension(2,2),intent(out):: vv,oo -!----------------------------------------------------------------------------- -real(dp):: a,b,c,d,e,f,g,h -!============================================================================= -a=em(1,1); b=em(1,2); c=em(2,2) -d=a*c-b*b! <- det(em) -e=(a+c)/2; f=(a-c)/2 -h=sqrt(f**2+b**2) -g=sqrt(b**2+(h+abs(f))**2) -if (g==0)then; vv(:,1)=(/u1,u0/) -elseif(f> 0)then; vv(:,1)=(/h+f,b/)/g -else ; vv(:,1)=(/b,h-f/)/g -endif -vv(:,2)=(/-vv(2,1),vv(1,1)/) -oo=matmul(transpose(vv),matmul(em,vv)) -oo(1,2)=u0; oo(2,1)=u0 -end subroutine eigensym2 -!============================================================================= -subroutine eigensym2d(em,vv,oo,vvd,ood,ff)! [eigensym2] -!============================================================================= -! For a symmetric 2*2 matrix, em, return the normalized eigenvectors, vv, and -! the diagonal matrix of eigenvalues, oo. If the two eigenvalues are equal, -! proceed no further and raise the logical failure flagg, ff, to .true.; -! otherwise, return with vvd=d(vv)/d(em) and ood=d(oo)/d(em) and ff=.false., -! and maintain the symmetries between the last two of the indices of -! these derivatives. -!============================================================================= -real(dp),dimension(2,2), intent(in ):: em -real(dp),dimension(2,2), intent(out):: vv,oo -real(dp),dimension(2,2,2,2),intent(out):: vvd,ood -logical, intent(out):: ff -!----------------------------------------------------------------------------- -real(dp),dimension(2,2):: emd,tt,vvr -real(dp) :: oodif,dtheta -integer :: i,j -!============================================================================= -call eigensym2(em,vv,oo); vvr(1,:)=-vv(2,:); vvr(2,:)=vv(1,:) -oodif=oo(1,1)-oo(2,2); ff=oodif==u0; if(ff)return -ood=0 -vvd=0 -do j=1,2 - do i=1,2 - emd=0 - if(i==j)then - emd(i,j)=u1 - else - emd(i,j)=o2; emd(j,i)=o2 - endif - tt=matmul(transpose(vv),matmul(emd,vv)) - ood(1,1,i,j)=tt(1,1) - ood(2,2,i,j)=tt(2,2) - dtheta=tt(1,2)/oodif - vvd(:,:,i,j)=vvr*dtheta - enddo -enddo -end subroutine eigensym2d - -!============================================================================= -subroutine invsym2(em,z)! [invsym2] -!============================================================================= -! Get the inverse of a 2*2 matrix (need not be symmetric in this case). -!============================================================================= -real(dp),dimension(2,2),intent(in ):: em -real(dp),dimension(2,2),intent(out):: z -!----------------------------------------------------------------------------- -real(dp):: detem -!============================================================================= -z(1,1)=em(2,2); z(2,1)=-em(2,1); z(1,2)=-em(1,2); z(2,2)=em(1,1) -detem=em(1,1)*em(2,2)-em(2,1)*em(1,2) -z=z/detem -end subroutine invsym2 -!============================================================================= -subroutine invsym2d(em,z,zd)! [invsym2] -!============================================================================= -! Get the inverse, z,of a 2*2 symmetric matrix, em, and its derivative, zd, -! with respect to symmetric variations of its components. I.e., for a -! symmetric infinitesimal change, delta_em, in em, the resulting -! infinitesimal change in z would be: -! delta_z(i,j) = matmul(zd(i,j,:,:),delta_em) -!============================================================================= -real(dp),dimension(2,2) ,intent(in ):: em -real(dp),dimension(2,2) ,intent(out):: z -real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- -integer:: k,l -!============================================================================= -call invsym2(em,z) -call id2222(zd) -do l=1,2; do k=1,2 - zd(:,:,k,l)=-matmul(matmul(z,zd(:,:,k,l)),z) -enddo; enddo -end subroutine invsym2d - -!============================================================================= -subroutine sqrtsym2(em,z)! [sqrtsym2] -!============================================================================= -! Get the sqrt of a symmetric positive-definite 2*2 matrix -!============================================================================= -real(dp),dimension(2,2),intent(in ):: em -real(dp),dimension(2,2),intent(out):: z -!----------------------------------------------------------------------------- -real(dp),dimension(2,2):: vv,oo -integer :: i -!============================================================================= -call eigensym2(em,vv,oo) -do i=1,2 -if(oo(i,i)<0)stop 'In sqrtsym2; matrix em is not non-negative' -oo(i,i)=sqrt(oo(i,i)); enddo -z=matmul(vv,matmul(oo,transpose(vv))) -end subroutine sqrtsym2 - -!============================================================================= -subroutine sqrtsym2d(x,z,zd)! [sqrtsym2] -!============================================================================= -! General routine to evaluate z=sqrt(x), and the symmetric -! derivative, zd = dz/dx, where x is a symmetric 2*2 positive-definite -! matrix. If the eigenvalues are very close together, extract their -! geometric mean for "preconditioning" a scaled version, px, of x, whose -! sqrt, and hence its derivative, can be easily obtained by the series -! expansion method. Otherwise, use the eigen-method (which entails dividing -! by the difference in the eignevalues to get zd, and which therefore -! fails when the eigenvalues become too similar). -!============================================================================= -real(dp),dimension(2,2), intent(in ):: x -real(dp),dimension(2,2), intent(out):: z -real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- -real(dp),dimension(2,2):: px -real(dp) :: rdetx,lrdetx,htrpxs,q -!============================================================================= -rdetx=sqrt(x(1,1)*x(2,2)-x(1,2)*x(2,1)) ! <- sqrt(determinant of x) -lrdetx=sqrt(rdetx) -px=x/rdetx ! <- preconditioned x (has unit determinant) -htrpxs= ((px(1,1)+px(2,2))/2)**2 ! <- {half-trace-px}-squared -q=htrpxs-u1 -if(q<.05)then ! <- Taylor expansion method - call sqrtsym2d_t(px,z,zd) - z=z*lrdetx; zd=zd/lrdetx -else - call sqrtsym2d_e(x,z,zd) ! <- Eigen-method -endif -end subroutine sqrtsym2d - -!============================================================================= -subroutine sqrtsym2d_e(x,z,zd)! [sqrtsym2d_e] -!============================================================================= -real(dp),dimension(2,2), intent(in ):: x -real(dp),dimension(2,2), intent(out):: z -real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- -real(dp),dimension(2,2,2,2):: vvd,ood -real(dp),dimension(2,2) :: vv,oo,oori,tt -integer :: i,j -logical :: ff -!============================================================================= -call eigensym2(x,vv,oo,vvd,ood,ff) -z=u0; z(1,1)=sqrt(oo(1,1)); z(2,2)=sqrt(oo(2,2)) -z=matmul(matmul(vv,z),transpose(vv)) -oori=0; oori(1,1)=u1/sqrt(oo(1,1)); oori(2,2)=u1/sqrt(oo(2,2)) -do j=1,2 -do i=1,2 - tt=matmul(vvd(:,:,i,j),transpose(vv)) - zd(:,:,i,j)=o2*matmul(matmul(matmul(vv,ood(:,:,i,j)),oori),transpose(vv))& - +matmul(tt,z)-matmul(z,tt) -enddo -enddo -end subroutine sqrtsym2d_e - -!============================================================================= -subroutine sqrtsym2d_t(x,z,zd)! [sqrtsym2d_t] -!============================================================================= -! Use the Taylor-series method (eigenvalues both fairly close to unity). -! For a 2*2 positive definite symmetric matrix x, try to get both the z=sqrt(x) -! and dz/dx using the binomial-expansion method applied to the intermediate -! matrix, r = (x-1). ie z=sqrt(x) = (1+r)^{1/2} = I + (1/2)*r -(1/8)*r^2 ... -! + [(-)^n *(2n)!/{(n+1)! * n! *2^{2*n-1}} ]*r^{n+1} -!============================================================================= -real(dp),dimension(2,2), intent(in ):: x -real(dp),dimension(2,2), intent(out):: z -real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- -integer,parameter :: nit=300 ! number of iterative increments allowed -real(dp),parameter :: crit=1.e-17 -real(dp),dimension(2,2) :: r,rp,rd,rpd -real(dp) :: c -integer :: i,j,n -!============================================================================= -r=x; r(1,1)=x(1,1)-1; r(2,2)=x(2,2)-1 -z=0; z(1,1)=u1; z(2,2)=u1 -rp=r -c=o2 -do n=0,nit - z=z+c*rp - rp=matmul(rp,r) - if(sum(abs(rp)) lx = ly = - a = - k = / diff --git a/ush/valid_param_vals.sh b/ush/valid_param_vals.sh index f0c78df05..d1eea61d5 100644 --- a/ush/valid_param_vals.sh +++ b/ush/valid_param_vals.sh @@ -47,8 +47,8 @@ valid_vals_CCPP_PHYS_SUITE=( \ ) valid_vals_OZONE_PARAM_NO_CCPP=("ozphys_2015" "ozphys") valid_vals_GFDLgrid_RES=("48" "96" "192" "384" "768" "1152" "3072") -valid_vals_EXTRN_MDL_NAME_ICS=("GSMGFS" "FV3GFS" "RAPX" "HRRRX") -valid_vals_EXTRN_MDL_NAME_LBCS=("GSMGFS" "FV3GFS" "RAPX" "HRRRX") +valid_vals_EXTRN_MDL_NAME_ICS=("GSMGFS" "FV3GFS" "RAPX" "HRRRX" "NAM") +valid_vals_EXTRN_MDL_NAME_LBCS=("GSMGFS" "FV3GFS" "RAPX" "HRRRX" "NAM") valid_vals_FV3GFS_FILE_FMT_ICS=("nemsio" "grib2") valid_vals_FV3GFS_FILE_FMT_LBCS=("nemsio" "grib2") valid_vals_GRID_GEN_METHOD=("GFDLgrid" "ESGgrid")