diff --git a/NUMERICS_PORT_VALIDATION.md b/NUMERICS_PORT_VALIDATION.md new file mode 100644 index 000000000..0309aef77 --- /dev/null +++ b/NUMERICS_PORT_VALIDATION.md @@ -0,0 +1,180 @@ +# Numerics/build port validation: `uob_fftw_sit23_numerics_2026` + +This documents validation of the FFTW3 spectral-transform core and HPC +build/environment support ported from `uob_fftw_sit23_new_maths2intel` onto +current `master`, done on branch `uob_fftw_sit23_numerics_2026`. Companion to +the physics port, validated separately in `PHYSICS_PORT_VALIDATION.md` on +`uob_fftw_sit23_physics_2026`. Two checks are covered: + +1. Does the standard (non-FFTW) build still produce bit-identical results to + `master`, given this branch touches several shared-physics files outside + any `#ifdef FFTW3` guard? +2. Does the new FFTW3 spectral core actually work — does it compile, link, + run, and produce physically correct output? + +## What was ported + +Branch `uob_fftw_sit23_numerics_2026` (commit `85bd52f1`) carries 9 commits +from `uob_fftw_sit23_new_maths2intel` onto current `master`, with original +authorship preserved (cherry-pick / diff+3-way-apply as appropriate): + +- **FFTW3 spectral transform core** (`grid_fourier.F90`, + `spherical_fourier.F90`, `spherical_fourier_fftw.F90` (new), `fft.F90`, + `fftw.F90` (new), plus `spectral_dynamics.F90`'s dispatch to it) — + George Lancaster, 2019. +- **`cb.enable_fftw3()` codebase API**, `DryCodeBaseFFTW` export, T213 + resolution, slurm multi-node support — sit23, 2020-2021. +- **`dry_fftw`/`grey_fftw` model directories** and the `held_suarez_fftw` + test case — qv18258, 2019. +- **Isambard/BluePebble/Cray-GNU HPC support**: env files, mkmf templates, + job scripts (kept per Stephen's explicit instruction — "just in case + they're useful for others" — even though these specific HPC systems may + be retired) and a set of compiler-compatibility fixes to shared physics + code needed to build on that toolchain — George Lancaster / Lancasterg, + 2019. +- **Postprocessing script improvements** — `run_plevel.py`'s `all_vars` + variable-subsetting option (reconciled against master's own since-modernised + scratch config), a `plevel_fn.py` fix for variables with a `scalar_axis` + dimension, and small fixes to `mppnccombine_run.sh` / + `compile_plev_interpolation.sh` / `modified_time_script.py` — sit23, + 2020-2021. + +**Deliberately excluded**, per Stephen's explicit decision: +- The branch's own incomplete single-precision attempt — deferred to the + separate, more complete `single_prec` branch, to be merged later as its + own project. +- `maths2` (env file + mkmf template) — turns out to **already be on + `master`**, added at commit `9e78f165` (this branch's own merge-base), so + it predates the branch entirely and was never part of the port. This + resolves the open question of whether `maths2` is the same lineage as the + currently-used `maths5`/`maths-gpu` nodes: it's moot, since master already + carries whatever `maths2` support it has independently of this branch. + +## The concern: two files change outside any `#ifdef FFTW3` guard + +Unlike the physics port (fully additive, all gated behind new namelist +options), this branch modifies shared code that runs on *every* build, +FFTW3 or not: + +- **`transforms.F90`**: `trans_spherical_to_grid_2d`/`_3d` and + `reverse_transpose_fourier` gain an optional `block` argument (defaults + preserve old behaviour); the `#ifdef INTERNAL_FILE_NML` namelist-reading + branch was rewritten to open `input.nml` directly via + `open_namelist_file()`/`close_file()` instead of reading from the + preloaded `input_nml_file` internal-file array. +- **`spectral_dynamics.F90`**: the same `INTERNAL_FILE_NML` namelist-reading + rewrite. +- The same `open_namelist_file()`/`close_file()` rewrite, applied + consistently, also appears in `diffusivity.F90`, `hs_forcing.F90`, + `lscale_cond.F90`, `monin_obukhov_kernel.F90`, `vert_turb_driver.F90`, + `vert_advection.F90`, `atmos_model.F90`, `atmosphere.F90`, + `polvani_2007.F90`, `spectral_initialize_fields.F90`, + `vert_coordinate.F90`, and `surface_flux.F90` (part of the Cray/Isambard + compiler-compatibility commit) — evidently a deliberate, systematic fix + for that toolchain, not an accidental one-off. + +None of these files add or rename any namelist variables (checked directly: +no diff in either module's `namelist /..._nml/` variable list), so unlike +the physics port's `update_land_mask_from_ice` case, `trip_test`'s +shared-namelist limitation ([[infra-trip-test-shared-namelist-limitation]]) +doesn't apply here — a direct `trip_test` comparison is valid. + +## Bit-reproducibility of the standard build + +Ran `trip_test_command_line` comparing `master` (`d1321ac3`) against this +branch's tip (`85bd52f1`) on `held_suarez` and `frierson` — chosen because +between them they exercise the dynamical core (`transforms.F90`, +`spectral_dynamics.F90`), the boundary-layer/turbulence physics +(`monin_obukhov_kernel.F90`, `vert_turb_driver.F90`, `diffusivity.F90`), +large-scale condensation (`lscale_cond.F90`), and the driver/coupler layer +(`atmos_model.F90`, `atmosphere.F90`, `surface_flux.F90`) — i.e. essentially +every non-FFTW-gated file this branch touches. + +**Result: pass on both.** + +``` +Test passed for held_suarez. Commit 85bd52f1 gives the same answer as commit d1321ac3 +Test passed for frierson. Commit 85bd52f1 gives the same answer as commit d1321ac3 +held_suarez : pass +frierson : pass +Congratulations, all tests have passed +``` + +Bit-identical output confirms the `open_namelist_file()` rewrite and the new +optional `block` argument are behaviourally neutral for the standard build — +consistent with what they look like on inspection (reading the same +`input.nml` content by a different mechanism; an optional argument that +defaults to the old behaviour everywhere it isn't explicitly passed). + +## Does the FFTW3 spectral core actually work? + +Compiling `held_suarez_fftw` with `-DFFTW3` was already confirmed to compile +cleanly during the port itself (both without and with FFTW3 enabled, the +latter needing the `dry_fftw`/`grey_fftw` `path_names` fix documented in the +commit history — stale manifests missing `cloud_simple`/ +`frierson_monin_obukhov.F90`). That only proves it compiles, not that it +runs correctly, so a further check was run here: **an actual 5-day +`held_suarez` integration**, comparing the FFTW3 spectral core against the +default (Temperton) FFT, same resolution (T42L25), same namelist, same +initial conditions. + +(This machine's `isca_env_26` conda environment doesn't ship FFTW3 dev +headers/libraries, unlike the target Isambard/BluePebble/Cray systems where +this would come from an environment module. A local-only, uncommitted +`-I`/`-L`/`-Wl,-rpath`/`-lfftw3` addition to `mkmf.template.ubuntu_conda`, +pointing at an existing FFTW3 install in the `isca_analysis` conda +environment, was used to make this testable here; reverted immediately +after, confirmed via `git status`/`git diff` afterward.) + +**Result: agrees with the reference FFT to floating-point roundoff, not +"close" — essentially exact.** + +| Field | max|diff| (FFTW3 − default FFT), day 5 | +|---|---| +| ps | 0.0 | +| temp | 0.0 | +| ucomp | 5.96×10-8 m/s | +| vcomp | 2.98×10-8 m/s | +| vor | 1.14×10-13 s-1 | + +These differences are all at or below double-precision rounding noise +(compare to typical field magnitudes: `ucomp`/`vcomp` range ±2–4 +m/s, `vor` range ±3×10-6 s-1) — several +orders of magnitude below the day-1 chaotic-noise-floor differences seen +between two independently-compiled binaries in the physics port's validation +(order 10-2 in physical units). Two different FFT algorithms +computing the same spherical harmonic transform, to this level of agreement +after 5 days of nonlinear integration, is strong evidence the FFTW3 dispatch +in `fft.F90`/`fftw.F90` is correctly wired up — not just "doesn't crash", +but numerically equivalent to the reference implementation it replaces. + +## Bottom line + +- Every file this branch touches outside an `#ifdef FFTW3` guard is + confirmed bit-reproducibility-neutral via `trip_test` (`held_suarez`, + `frierson`, both pass against `master`). +- The FFTW3 spectral core, once actually run (not just compiled), agrees + with the reference FFT implementation to floating-point roundoff over a + 5-day integration — the transform is correctly implemented, not just + syntactically valid. +- Both the standard and FFTW3-enabled builds compile cleanly, including + after fixing a real bug surfaced during the port itself (stale + `dry_fftw`/`grey_fftw` `path_names` manifests missing modules current + master's `dry`/`grey` models depend on). +- `maths2` and single-precision support are out of scope for this branch + (see above) — nothing further needed here for either. + +## Appendix: how this was tested + +`trip_test_command_line -e held_suarez frierson -n 8 -r + d1321ac3 85bd52f1`, run locally against this +worktree as the `-r` repo argument (a plain local git clone source — no +GitHub round-trip needed, since all worktrees on this machine share one +`.git` object store). + +The FFTW3-vs-default-FFT comparison used a standalone script (not +`trip_test`, since it needs to compile two different `CodeBase` subclasses — +`DryCodeBaseFFTW` vs `DryCodeBase` — against the *same* commit, which +`trip_test` isn't set up to do): both compiled from this branch's worktree, +run for 5 days at T42L25 with identical namelists/diagnostics/initial +conditions, output compared field-by-field with `xarray`. diff --git a/exp/site_specific/bristol_bcp3/README.md b/exp/site_specific/bristol_bcp3/README.md new file mode 100644 index 000000000..c1dfac65b --- /dev/null +++ b/exp/site_specific/bristol_bcp3/README.md @@ -0,0 +1,4 @@ +# Instructions for running Isca on BlueCrystal 3 at Bristol + +These instructions are intended to get you up-and-running with a simple Held-Suarez test case. They assume you are starting with a default user environment on BlueCrystal 3, so some changes might be needed if you have already modified your environment. + diff --git a/exp/site_specific/bristol_bcp4/ReadMe_Bristol.md b/exp/site_specific/bristol_bcp4/ReadMe_Bristol.md new file mode 100644 index 000000000..5928065c4 --- /dev/null +++ b/exp/site_specific/bristol_bcp4/ReadMe_Bristol.md @@ -0,0 +1,73 @@ + +# Instructions for running Isca on BlueCrystal 4 at Bristol + +These instructions are intended to get you up-and-running with a simple Held-Suarez test case. They assume you are starting with a default user environment on BlueCrystal 4, so some changes might be needed if you have already modified your environment. Please don't hesistate to get in touch (w.seviour@bristol.ac.uk) if you have any questions. + +First, to enable `git`, run: + +```{bash} +$ module load git +``` + +Now you can clone the Isca repository: + +```{bash} +$ git clone git@github.com:ExeClim/Isca.git +$ cd Isca +``` + +Before you can run Isca, you'll need to load the Anaconda python distribution: + +```{bash} +$ module load languages/anaconda3 +``` + +Next we'll make a conda environment for Isca (this means it will have all the right versions of the various packages on which it depends): + +```{bash} +$ conda create -n isca_env python ipython +$ source activate isca_env +(isca_env) $ cd Isca/src/extra/python +(isca_env) $ pip install -r requirements.txt + +Successfully installed MarkupSafe-1.0 f90nml jinja2-2.9.6 numpy-1.13.3 pandas-0.21.0 python-dateutil-2.6.1 pytz-2017.3 sh-1.12.14 six-1.11.0 xarray-0.9.6 + +(isca_env) $ pip install -e . +... +Successfully installed Isca +``` + +Finally, we'll need to update the `~/.bashrc` file. Add the following lines: + +```{bash} +# directory of the Isca source code +export GFDL_BASE=$HOME/Isca +# "environment" configuration for bc4 +export GFDL_ENV=bristol-bc4 +# temporary working directory used in running the model +export GFDL_WORK=/mnt/storage/home/$USER/scratch/Isca_work +# directory for storing model output +export GFDL_DATA=/mnt/storage/home/$USER/scratch/Isca_data +``` + +Then make the `Isca_work` and `Isca_data` directories: + +```{bash} +(isca_env) $ mkdir -p /mnt/storage/home/$USER/scratch/Isca_work +(isca_env} $ mkdir -p /mnt/storage/home/$USER/scratch/Isca_data +(isca_env) $ bash +``` +Now everything should be set up and we can try a test run. The following should compile and run 12 months of a Held-Suarez test case, at T42 resolution spread over 16 cores. + +```{bash} +(isca_env) $ cd Isca/exp/site-specific +(isca_env) $ sbatch isca_slurm_job.sh +``` + +This should produce a slurm file showing the progress as it compiles and runs. To track the progress: + +```{bash} +(isca_env) $ tail -f slurm_*.o +``` + +All being well, after about 20 minutes the job should complete, and you'll find some output files in `/mnt/storage/home/$USER/scratch/Isca_data`. \ No newline at end of file diff --git a/exp/site_specific/bristol_bcp4/isca_slurm_job.sh b/exp/site_specific/bristol_bcp4/isca_slurm_job.sh new file mode 100644 index 000000000..d3850f522 --- /dev/null +++ b/exp/site_specific/bristol_bcp4/isca_slurm_job.sh @@ -0,0 +1,25 @@ +#!/bin/bash -l + +#SBATCH --job-name=held_suarez_test_case +#SBATCH --partition=veryshort +#SBATCH --time=1:00:00 +#SBATCH --nodes=1 +#number of tasks ~ processes per node +#SBATCH --ntasks-per-node=16 +#number of cpus (cores) per task (process) +#SBATCH --cpus-per-task=1 +#SBATCH --output=slurm_%j.o + +echo Running on host `hostname` +echo Time is `date` +echo Directory is `pwd` + + +module purge +source $HOME/.bashrc +source $GFDL_BASE/src/extra/env/bristol-bc4 +source activate isca_env + + + +$HOME/.conda/envs/isca_env/bin/python $GFDL_BASE/exp/test_cases/held_suarez/held_suarez_test_case.py diff --git a/exp/site_specific/bristol_bluepebble/README.md b/exp/site_specific/bristol_bluepebble/README.md new file mode 100644 index 000000000..81ccbad2c --- /dev/null +++ b/exp/site_specific/bristol_bluepebble/README.md @@ -0,0 +1,73 @@ +# Instructions for running Isca on Bristol's BluePebble Supercomputer (Intel compiler) + +These instructions are intended to get you up-and-running with a simple Held-Suarez test case. They assume you are starting with a default user environment on BluePebble, so some changes might be needed if you have already modified your environment. + +First, you need to clone the Isca repository into your home directory + +```{bash} +$ git clone git@github.com:ExeClim/Isca.git +$ cd Isca +``` + +Before you can run Isca, you'll need to load the CrayPython module: + +```{bash} +$ module load cray-python/3.6.5.6 +``` + +Next we'll make a Python environment for Isca (this means it will have all the right versions of the various packages on which it depends): + +```{bash} +$ python3 -m venv ~/isca_env +$ source ~/isca_env/bin/activate +(isca_env) $ cd Isca/src/extra/python +(isca_env) $ pip install -r requirements.txt + +Successfully installed MarkupSafe-1.0 f90nml jinja2-2.9.6 numpy-1.13.3 pandas-0.21.0 python-dateutil-2.6.1 pytz-2017.3 sh-1.12.14 six-1.11.0 xarray-0.9.6 + +(isca_env) $ pip install -e . +... +Successfully installed Isca +``` + +Finally, we'll need to update the `~/.bash_profile` file. Add the following lines: + +```{bash} +# directory of the Isca source code +export GFDL_BASE=$HOME/Isca +# "environment" configuration for bc4 +export GFDL_ENV=bluepebble +# temporary working directory used in running the model +export GFDL_WORK=$HOME/Isca_work +# directory for storing model output +export GFDL_DATA=$HOME/Isca_data +``` + +Then make the `Isca_work` and `Isca_data` directories: + +```{bash} +(isca_env) $ mkdir -p $HOME/isca_home/isca_work +(isca_env} $ mkdir -p $HOME/isca_home/isca_data +(isca_env) $ bash +``` + +Now everything should be set up and we can try a test run. The following should compile and run 12 months of a Held-Suarez test case. + +```{bash} +(isca_env) $ cd Isca/exp/site-specific/bristol_bluepebble +(isca_env) $ qsub grey +``` + +This should produce a two outputfiles: bluepebble.job.e and bluepebble.job.e reporting the standard output and error output of the job. + +You can view details about the job using + +```{bash} +(isca_env) $ qstat -f +``` + +All being well, after about 60 minutes the job should complete, and you'll find some output files in `$GFDL_DATA`. + +```{bash} +(isca_env) $ cd $GFDL_DATA +``` diff --git a/exp/site_specific/bristol_bluepebble/bluepebble.job b/exp/site_specific/bristol_bluepebble/bluepebble.job new file mode 100644 index 000000000..df423e5bb --- /dev/null +++ b/exp/site_specific/bristol_bluepebble/bluepebble.job @@ -0,0 +1,8 @@ +#!/bin/sh +#PBS -l select=1:ncpus=16:mem=20GB +#PBS -l walltime=04:00:00 + +module load tools/git/2.22.0 +source activate isca_env + +python $GFDL_BASE/exp/test_cases/held_suarez/held_suarez_test_case.py diff --git a/exp/site_specific/gw4_isambard/README.md b/exp/site_specific/gw4_isambard/README.md new file mode 100644 index 000000000..6ab04e634 --- /dev/null +++ b/exp/site_specific/gw4_isambard/README.md @@ -0,0 +1,80 @@ +# Instructions for running Isca on GW4's Isambard Supercomputer (GNU compiler) + +These instructions are intended to get you up-and-running with a simple Held-Suarez test case. They assume you are starting with a default user environment on Isambard, so some changes might be needed if you have already modified your environment. + +Isamabard is a cray supercomputer comprised of Arm processor nodes. This means that configuration is slightly different to that of the usual Intel-based machines. + +First, you need to clone the Isca repository into your home directory + +```{bash} +$ git clone git@github.com:ExeClim/Isca.git +$ cd Isca +``` + +Before you can run Isca, you'll need to load the CrayPython module: + +```{bash} +$ module load cray-python/3.6.5.6 +``` + +Next we'll make a Python environment for Isca (this means it will have all the right versions of the various packages on which it depends): + +```{bash} +$ python3 -m venv ~/isca_env +$ source ~/isca_env/bin/activate +(isca_env) $ cd Isca/src/extra/python +(isca_env) $ pip install -r requirements.txt + +Successfully installed MarkupSafe-1.0 f90nml jinja2-2.9.6 numpy-1.13.3 pandas-0.21.0 python-dateutil-2.6.1 pytz-2017.3 sh-1.12.14 six-1.11.0 xarray-0.9.6 + +(isca_env) $ pip install -e . +... +Successfully installed Isca +``` + +Finally, we'll need to update the `~/.bashrc` file. Add the following lines: + +```{bash} +# directory of the Isca source code +export GFDL_BASE=$HOME/Isca +# "environment" configuration for bc4 +export GFDL_ENV=bristol-bc4 +# temporary working directory used in running the model +export GFDL_WORK=$HOME/Isca_work +# directory for storing model output +export GFDL_DATA=$HOME/Isca_data + +#Cray-specific environment variables +export EXECUTION_TYPE=APRUN +#export GFDL_MKMF_TEMPLATE=isambard_gfortran #export GFDL_ENV=isambard_gfortran +module unload PrgEnv-cray/6.0.5 +``` + +Then make the `Isca_work` and `Isca_data` directories: + +```{bash} +(isca_env) $ mkdir -p $HOME/isca_home/isca_work +(isca_env} $ mkdir -p $HOME/isca_home/isca_data +(isca_env) $ bash +``` + +Now everything should be set up and we can try a test run. The following should compile and run 12 months of a Held-Suarez test case, at T85 resolution spread over 64 cores. + +```{bash} +(isca_env) $ cd Isca/exp/site-specific/gw4_isambard +(isca_env) $ qsub isca_pbspro.job +``` + +This should produce a two outputfiles: isca_pbspro.job.e and isca_pbspro.job.e reporting the standard output and error output of the job. + +You can view details about the job using + +```{bash} +(isca_env) $ qstat -f +``` + +All being well, after about 60 minutes the job should complete, and you'll find some output files in `$GFDL_DATA`. + +```{bash} +(isca_env) $ cd $GFDL_DATA +``` diff --git a/exp/site_specific/gw4_isambard/held_suarez_test_case.py b/exp/site_specific/gw4_isambard/held_suarez_test_case.py new file mode 100644 index 000000000..fe426a706 --- /dev/null +++ b/exp/site_specific/gw4_isambard/held_suarez_test_case.py @@ -0,0 +1,109 @@ +import numpy as np + +from isca import DryCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE + +NCORES = 64 +RESOLUTION = 'T85', 25 # T42 horizontal resolution, 25 levels in pressure + +# a CodeBase can be a directory on the computer, +# useful for iterative development +cb = DryCodeBase.from_directory(GFDL_BASE) + +# or it can point to a specific git repo and commit id. +# This method should ensure future, independent, reproducibility of results. +# cb = DryCodeBase.from_repo(repo='https://github.com/isca/isca', commit='isca1.1') + +# compilation depends on computer specific settings. The $GFDL_ENV +# environment variable is used to determine which `$GFDL_BASE/src/extra/env` file +# is used to load the correct compilers. The env file is always loaded from +# $GFDL_BASE and not the checked out git repo. + +cb.compile() # compile the source code to working directory $GFDL_WORK/codebase + +# create an Experiment object to handle the configuration of model parameters +# and output diagnostics + +exp_name = 'held_suarez_default' +exp = Experiment(exp_name, codebase=cb) + +#Tell model how to write diagnostics +diag = DiagTable() +diag.add_file('atmos_monthly', 30, 'days', time_units='days') + +#Tell model which diagnostics to write +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk') +diag.add_field('dynamics', 'pk') +diag.add_field('dynamics', 'ucomp', time_avg=True) +diag.add_field('dynamics', 'vcomp', time_avg=True) +diag.add_field('dynamics', 'temp', time_avg=True) +diag.add_field('dynamics', 'vor', time_avg=True) +diag.add_field('dynamics', 'div', time_avg=True) + +exp.diag_table = diag + +# define namelist values as python dictionary +# wrapped as a namelist object. +namelist = Namelist({ + 'main_nml': { + 'dt_atmos': 600, + 'days': 30, + 'calendar': 'thirty_day', + 'current_date': [2000,1,1,0,0,0] + }, + + 'atmosphere_nml': { + 'idealized_moist_model': False # False for Newtonian Cooling. True for Isca/Frierson + }, + + 'spectral_dynamics_nml': { + 'damping_order' : 4, # default: 2 + 'water_correction_limit' : 200.e2, # default: 0 + 'reference_sea_level_press': 1.0e5, # default: 101325 + 'valid_range_t' : [100., 800.], # default: (100, 500) + 'initial_sphum' : 0.0, # default: 0 + 'vert_coord_option' : 'uneven_sigma', # default: 'even_sigma' + 'scale_heights': 6.0, + 'exponent': 7.5, + 'surf_res': 0.5 + }, + + # configure the relaxation profile + 'hs_forcing_nml': { + 't_zero': 315., # temperature at reference pressure at equator (default 315K) + 't_strat': 200., # stratosphere temperature (default 200K) + 'delh': 60., # equator-pole temp gradient (default 60K) + 'delv': 10., # lapse rate (default 10K) + 'eps': 0., # stratospheric latitudinal variation (default 0K) + 'sigma_b': 0.7, # boundary layer friction height (default p/ps = sigma = 0.7) + + # negative sign is a flag indicating that the units are days + 'ka': -40., # Constant Newtonian cooling timescale (default 40 days) + 'ks': -4., # Boundary layer dependent cooling timescale (default 4 days) + 'kf': -1., # BL momentum frictional timescale (default 1 days) + + 'do_conserve_energy': True, # convert dissipated momentum into heat (default True) + }, + + 'diag_manager_nml': { + 'mix_snapshot_average_fields': False + }, + + 'fms_nml': { + 'domains_stack_size': 600000 # default: 0 + }, + + 'fms_io_nml': { + 'threading_write': 'single', # default: multi + 'fileset_write': 'single', # default: multi + } +}) + +exp.namelist = namelist +exp.set_resolution(*RESOLUTION) + +#Lets do a run! +if __name__ == '__main__': + exp.run(1, num_cores=NCORES, use_restart=False) + for i in range(2, 13): + exp.run(i, num_cores=NCORES) # use the restart i-1 by default \ No newline at end of file diff --git a/exp/site_specific/gw4_isambard/isca_pbspro.job b/exp/site_specific/gw4_isambard/isca_pbspro.job new file mode 100644 index 000000000..d432bf2f7 --- /dev/null +++ b/exp/site_specific/gw4_isambard/isca_pbspro.job @@ -0,0 +1,8 @@ +#!/bin/sh +#PBS -q arm +#PBS -l select=1:ncpus=64 +#PBS -l walltime=01:00:00 + +source ~/pyi3/bin/activate +python $GFDL_BASE/exp/site_specific/gw4_isambard/held_suarez_test_case.py + diff --git a/exp/test_cases/held_suarez_fftw/held_suarez_test_case_fftw.py b/exp/test_cases/held_suarez_fftw/held_suarez_test_case_fftw.py new file mode 100644 index 000000000..0d1903fb8 --- /dev/null +++ b/exp/test_cases/held_suarez_fftw/held_suarez_test_case_fftw.py @@ -0,0 +1,111 @@ +import numpy as np + +from isca import DryCodeBaseFFTW, DiagTable, Experiment, Namelist, GFDL_BASE + +NCORES = 2 +RESOLUTION = 'T42', 25 # T42 horizontal resolution, 25 levels in pressure + +# a CodeBase can be a directory on the computer, +# useful for iterative development +cb = DryCodeBaseFFTW.from_directory(GFDL_BASE) + +# or it can point to a specific git repo and commit id. +# This method should ensure future, independent, reproducibility of results. +# cb = DryCodeBase.from_repo(repo='https://github.com/isca/isca', commit='isca1.1') + +# compilation depends on computer specific settings. The $GFDL_ENV +# environment variable is used to determine which `$GFDL_BASE/src/extra/env` file +# is used to load the correct compilers. The env file is always loaded from +# $GFDL_BASE and not the checked out git repo. + +# compile using fftw3 +cb.compile_flags.append('-DFFTW3') +cb.compile() # compile the source code to working directory $GFDL_WORK/codebase + +# create an Experiment object to handle the configuration of model parameters +# and output diagnostics + +exp_name = 'held_suarez_default' +exp = Experiment(exp_name, codebase=cb) + +#Tell model how to write diagnostics +diag = DiagTable() +diag.add_file('atmos_monthly', 30, 'days', time_units='days') + +#Tell model which diagnostics to write +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk') +diag.add_field('dynamics', 'pk') +diag.add_field('dynamics', 'ucomp', time_avg=True) +diag.add_field('dynamics', 'vcomp', time_avg=True) +diag.add_field('dynamics', 'temp', time_avg=True) +diag.add_field('dynamics', 'vor', time_avg=True) +diag.add_field('dynamics', 'div', time_avg=True) + +exp.diag_table = diag + +# define namelist values as python dictionary +# wrapped as a namelist object. +namelist = Namelist({ + 'main_nml': { + 'dt_atmos': 600, + 'days': 30, + 'calendar': 'thirty_day', + 'current_date': [2000,1,1,0,0,0] + }, + + 'atmosphere_nml': { + 'idealized_moist_model': False # False for Newtonian Cooling. True for Isca/Frierson + }, + + 'spectral_dynamics_nml': { + 'damping_order' : 4, # default: 2 + 'water_correction_limit' : 200.e2, # default: 0 + 'reference_sea_level_press': 1.0e5, # default: 101325 + 'valid_range_t' : [100., 800.], # default: (100, 500) + 'initial_sphum' : 0.0, # default: 0 + 'vert_coord_option' : 'uneven_sigma', # default: 'even_sigma' + 'scale_heights': 6.0, + 'exponent': 7.5, + 'surf_res': 0.5 + }, + + # configure the relaxation profile + 'hs_forcing_nml': { + 't_zero': 315., # temperature at reference pressure at equator (default 315K) + 't_strat': 200., # stratosphere temperature (default 200K) + 'delh': 60., # equator-pole temp gradient (default 60K) + 'delv': 10., # lapse rate (default 10K) + 'eps': 0., # stratospheric latitudinal variation (default 0K) + 'sigma_b': 0.7, # boundary layer friction height (default p/ps = sigma = 0.7) + + # negative sign is a flag indicating that the units are days + 'ka': -40., # Constant Newtonian cooling timescale (default 40 days) + 'ks': -4., # Boundary layer dependent cooling timescale (default 4 days) + 'kf': -1., # BL momentum frictional timescale (default 1 days) + + 'do_conserve_energy': True, # convert dissipated momentum into heat (default True) + }, + + 'diag_manager_nml': { + 'mix_snapshot_average_fields': False + }, + + 'fms_nml': { + 'domains_stack_size': 600000 # default: 0 + }, + + 'fms_io_nml': { + 'threading_write': 'single', # default: multi + 'fileset_write': 'single', # default: multi + } +}) + +exp.namelist = namelist +exp.set_resolution(*RESOLUTION) + +#Lets do a run! +if __name__ == '__main__': + exp.run(1, num_cores=NCORES, use_restart=False) + for i in range(2, 13): + exp.run(i, num_cores=NCORES) # use the restart i-1 by default diff --git a/exp/test_cases/held_suarez_fftw/readme.md b/exp/test_cases/held_suarez_fftw/readme.md new file mode 100644 index 000000000..c662bcb2c --- /dev/null +++ b/exp/test_cases/held_suarez_fftw/readme.md @@ -0,0 +1,37 @@ +# Held-Suarez using FFTW + +This configuration is idential to the Held-Suarez configuration found in `../held_suarez`, however it uses the FFTW library to perform the FFT instead of the default Temperton's FFT. + +## Running this configuration +To run this configuration, you need to ensure that you have FFTW installed on your system, and that it can be linked to the Isca code. This usually involves loading the relevant module file. For BCP4 at Bristol, this is done by adding the following command to the environment file: + +``` +module load libs/fftw/3.3.6 +``` + +Once this has been done, you can simply run from within this folder: + +``` +python held_suarez_test_case_fftw.py +``` + +## Other configurations +In order to use FFTW with other model configurations there are two steps you must take: + +### 1. Add the correct preprocessor directive + +You must add the `-DFFTW3` preprocessor directive to the list of compiler flags that will be used to compile Isca. This can be done by adding the following line to your python script: + +``` +cb.compile_flags.append('-DFFTW3') +``` +This must be done **before** calling `cb.compile()`. An example of this is given in the `held_suarez_test_case_fftw.py` file inside this folder. + +### 2. Add the fftw.F90 module + +To use FFTW, a wrapper module has been created (`Isca/src/shared/fft/fftw.F90`). This file must be added to the list of files that are to be compiled using your configuration, as done in `src/extra/model/isca_fftw/path_names`. + + + +Having followed both steps, the FFTW3 library should be configured, however feel free to send me an email at qv18258@bristol.ac.uk if you get stuck. + diff --git a/exp/test_cases/isca_job_bluepebble.job b/exp/test_cases/isca_job_bluepebble.job new file mode 100644 index 000000000..303179131 --- /dev/null +++ b/exp/test_cases/isca_job_bluepebble.job @@ -0,0 +1,8 @@ +#!/bin/sh +#PBS -q mvlong +#PBS -l select=1:ncpus=36 +#PBS -l walltime=0:60:00 + +module load lang/python/anaconda/3.7-2019.03.biopython +source activate isca_env +python held_suarez/held_suarez_test_case.py diff --git a/postprocessing/mppnccombine_run.sh b/postprocessing/mppnccombine_run.sh index d43e7bafc..792858e7e 100755 --- a/postprocessing/mppnccombine_run.sh +++ b/postprocessing/mppnccombine_run.sh @@ -1,10 +1,13 @@ #!/usr/bin/env bash -# 1. Load the necessary tools into the environment -module purge +# 1. If safe to purge, then do so +if [[ -z "${GFDL_PURGE}" ]] || [ "${GFDL_PURGE}" = true ]; then + module purge +fi +# Load the necessary tools into the environment source $GFDL_BASE/src/extra/env/$GFDL_ENV module list # Run with input file as input taken from BASH script echo ${1} -exec ${1}/mppnccombine.x ${2} \ No newline at end of file +exec ${1}/mppnccombine.x ${2} diff --git a/postprocessing/plevel_interpolation/compile_plev_interpolation.sh b/postprocessing/plevel_interpolation/compile_plev_interpolation.sh index c9c668813..fa2790dde 100755 --- a/postprocessing/plevel_interpolation/compile_plev_interpolation.sh +++ b/postprocessing/plevel_interpolation/compile_plev_interpolation.sh @@ -5,6 +5,9 @@ cd ./exec source $GFDL_BASE/src/extra/env/$GFDL_ENV -../bin/mkmf -p plev.x -t ../bin/mkmf.template.ia64 -c "-Duse_netCDF" -a ../src ../src/path_names ../src/shared/mpp/include ../src/shared/include +compiler=${GFDL_MKMF_TEMPLATE:-ia64} +template=mkmf.template.${compiler} + +../bin/mkmf -p plev.x -t $GFDL_BASE/src/extra/python/isca/templates/$template -c "-Duse_netCDF" -a ../src ../src/path_names ../src/shared/mpp/include ../src/shared/include make -f Makefile diff --git a/postprocessing/plevel_interpolation/scripts/plevel_fn.py b/postprocessing/plevel_interpolation/scripts/plevel_fn.py index 2b7ae3f02..b6b70423a 100644 --- a/postprocessing/plevel_interpolation/scripts/plevel_fn.py +++ b/postprocessing/plevel_interpolation/scripts/plevel_fn.py @@ -80,7 +80,6 @@ def two_daily_average(nc_file_in, nc_file_out, avg_or_daily): subprocess.call('cdo timselmean,'+str(number_of_timesteps)+' '+nc_file_in+' '+nc_file_out, shell=True) def join_files(files_in, file_name_out): - subprocess.call('cdo mergetime '+files_in+' '+file_name_out, shell=True) def join_files_base_dir(base_dir, files_in, file_name_out): diff --git a/postprocessing/plevel_interpolation/scripts/run_plevel.py b/postprocessing/plevel_interpolation/scripts/run_plevel.py index b6a7cfce6..6059a4836 100644 --- a/postprocessing/plevel_interpolation/scripts/run_plevel.py +++ b/postprocessing/plevel_interpolation/scripts/run_plevel.py @@ -5,22 +5,27 @@ import time import pdb import subprocess +import numpy as np start_time=time.time() -# base_dir='/disca/share/sit204/data_from_isca_cpu/cssp_perturb_exps/anoms/' base_dir = os.environ['GFDL_DATA'] -#exp_name_list = ['soc_ga3_files_smooth_topo_fftw_mk1_fresh_compile_long', 'soc_ga3_files_smooth_topo_old_fft_mk2_long'] exp_name_list = ['ml_test_with_ml_full_sd_mean_clim_half_month'] avg_or_daily_list=['half_monthly'] start_file=133 end_file=134 +all_vars=True +var_names_list = 'slp height precipitation vcomp ucomp temp_2m temp div flux_t flux_lhe bucket_depth' nfiles=(end_file-start_file)+1 do_extra_averaging=False #If true, then 6hourly data is averaged into daily data using cdo -group_months_into_one_file=False # If true then monthly data files and daily data files are merged into one big netcdf file each. +group_months_into_one_file=True # If true then monthly data files and daily data files are merged into one big netcdf file each. +overwrite_previous_combined_files=True +n_splits_of_combined_nc_file = 20 level_set='standard' #Default is the standard levels used previously. ssw_diagnostics are the ones blanca requested for MiMa validation mask_below_surface_set=' ' #Default is to mask values that lie below the surface pressure when interpolated. For some applications, e.g. Tom Clemo's / Mark Baldwin's stratosphere index, you want to have values interpolated below ground, i.e. as if the ground wasn't there. To use this option, this value should be set to '-x '. - +all_vars=True +#var_names_list = 'slp height precipitation vcomp ucomp temp_2m temp div flux_t flux_lhe bucket_depth' +var_names_list = '-a slp height' try: out_dir @@ -43,13 +48,22 @@ plevs['6hourly']=' -p "1000 10000 25000 50000 85000 92500"' plevs['daily'] =' -p "1000 10000 25000 50000 85000 92500"' - var_names['monthly']='-a slp height' - var_names['half_monthly']='-a slp height' - var_names['pentad']='-a slp height' + if all_vars: + var_names['monthly']='-a slp height' + var_names['half_monthly']='-a slp height' + var_names['pentad']='-a slp height' + else: + var_names['monthly']=var_names_list + var_names['half_monthly']=var_names_list + var_names['pentad']=var_names_list var_names['timestep']='-a' var_names['6hourly']='ucomp slp height vor t_surf vcomp omega' var_names['daily']='ucomp slp height vor t_surf vcomp omega temp' - file_suffix='_interp_new_height_temp_not_below_ground' + + if all_vars: + file_suffix='_interp_new_height_temp_not_below_ground' + else: + file_suffix='_interp_new_height_temp_not_below_ground_subset_vars_5' elif level_set=='ssw_diagnostics': plevs['6hourly']=' -p "1000 10000"' @@ -85,18 +99,25 @@ # two_daily_average(nc_file_out, nc_file_out_two_daily, avg_or_daily) if group_months_into_one_file: - avg_or_daily_list_together=['daily'] + avg_or_daily_list_together=avg_or_daily_list for exp_name in exp_name_list: for avg_or_daily in avg_or_daily_list_together: - nc_file_string='' - for n in range(nfiles): - nc_file_in = base_dir+'/'+exp_name+'/run'+number_prefix+str(n+start_file)+'/atmos_'+avg_or_daily+file_suffix+'.nc' - nc_file_string=nc_file_string+' '+nc_file_in - nc_file_out=base_dir+'/'+exp_name+'/atmos_'+avg_or_daily+'_together'+file_suffix+'.nc' - if not os.path.isfile(nc_file_out): - join_files(nc_file_string,nc_file_out) + + for itr_idx in range(n_splits_of_combined_nc_file): + nc_file_string='' + n_files_in_each = int(np.ceil(nfiles / n_splits_of_combined_nc_file)) + + start_file_itr = start_file + (itr_idx)*n_files_in_each + end_file_itr = int(np.min([start_file + ((itr_idx+1)*n_files_in_each - 1), end_file])) + + for n in range(start_file_itr, end_file_itr+1): + nc_file_in = base_dir+'/'+exp_name+'/run%04d'%(n)+'/atmos_'+avg_or_daily+file_suffix+'.nc' + nc_file_string=nc_file_string+' '+nc_file_in + nc_file_out=base_dir+'/'+exp_name+'/atmos_'+avg_or_daily+f'_together_{start_file_itr}_{end_file_itr}'+file_suffix+'.nc' + if not os.path.isfile(nc_file_out) or overwrite_previous_combined_files: + join_files(nc_file_string,nc_file_out) print('execution time', time.time()-start_time) diff --git a/src/atmos_param/diffusivity/diffusivity.F90 b/src/atmos_param/diffusivity/diffusivity.F90 index 921e5bb73..f79aef19e 100644 --- a/src/atmos_param/diffusivity/diffusivity.F90 +++ b/src/atmos_param/diffusivity/diffusivity.F90 @@ -180,7 +180,9 @@ subroutine diffusivity_init if (file_exist('input.nml')) then #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=diffusivity_nml, iostat=io) + unit = open_namelist_file() + read (unit, nml=diffusivity_nml, iostat=io) + call close_file(unit) ierr = check_nml_error(io,"diffusivity_nml") #else unit = open_namelist_file () diff --git a/src/atmos_param/hs_forcing/hs_forcing.F90 b/src/atmos_param/hs_forcing/hs_forcing.F90 index 164e82f36..8dd7cb0a9 100644 --- a/src/atmos_param/hs_forcing/hs_forcing.F90 +++ b/src/atmos_param/hs_forcing/hs_forcing.F90 @@ -36,7 +36,8 @@ module hs_forcing_mod check_nml_error, & mpp_pe, mpp_root_pe, close_file, & write_version_number, stdlog, & - uppercase, read_data, write_data, set_domain + uppercase, read_data, write_data, & + set_domain, open_namelist_file use time_manager_mod, only: time_type, get_time @@ -306,7 +307,9 @@ subroutine hs_forcing_init ( axes, Time, lonb, latb, lat ) ! ----- read namelist ----- #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=hs_forcing_nml, iostat=io) + unit = open_namelist_file ( ) + read (unit, nml=hs_forcing_nml, iostat=io) + call close_file(unit) ierr = check_nml_error(io, 'hs_forcing_nml') #else if (file_exist('input.nml')) then diff --git a/src/atmos_param/lscale_cond/lscale_cond.F90 b/src/atmos_param/lscale_cond/lscale_cond.F90 index 231c03803..14a38c779 100644 --- a/src/atmos_param/lscale_cond/lscale_cond.F90 +++ b/src/atmos_param/lscale_cond/lscale_cond.F90 @@ -264,13 +264,15 @@ subroutine lscale_cond_init () ! !----------------------------------------------------------------------- - integer unit,io,ierr, logunit + integer unit, io, ierr, logunit !----------- read namelist --------------------------------------------- #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=lscale_cond_nml, iostat=io) - ierr = check_nml_error(io,"lscale_cond_nml") + unit = open_namelist_file ( ) + read (unit, nml=lscale_cond_nml, iostat=io) + call close_file(unit) + ierr = check_nml_error(io,"lscale_cond_nml") #else if (file_exist('input.nml')) then unit = open_namelist_file () diff --git a/src/atmos_param/monin_obukhov/monin_obukhov_kernel.F90 b/src/atmos_param/monin_obukhov/monin_obukhov_kernel.F90 index 45b9db986..2e95d5e45 100644 --- a/src/atmos_param/monin_obukhov/monin_obukhov_kernel.F90 +++ b/src/atmos_param/monin_obukhov/monin_obukhov_kernel.F90 @@ -19,6 +19,37 @@ !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! User must define n based on model resolution +!| Resolution | DIM | +!|------------|-----| +!| T21 | 64 | +!| T42 | 128 | +!| T85 | 256 | +!| T170 | 512 | + +#define CRAY_DIM n + +#ifdef cray_dim_64 +#undef CRAY_DIM +#define CRAY_DIM 64 +#endif + +#ifdef cray_dim_128 +#undef CRAY_DIM +#define CRAY_DIM 128 +#endif + +#ifdef cray_dim_256 +#undef CRAY_DIM +#define CRAY_DIM 256 +#endif + +#ifdef cray_dim_512 +#undef CRAY_DIM +#define CRAY_DIM 512 +#endif + + #include ! -*-F90-*- @@ -143,7 +174,7 @@ _PURE subroutine monin_obukhov_drag_1d(grav, vonkarm, & logical, intent(in ), dimension(n) :: avail ! provided mask integer, intent(out ) :: ier - real , dimension(n) :: rich, fm, ft, fq, zz + real , dimension(CRAY_DIM) :: rich, zz, fm, ft, fq logical, dimension(n) :: mask, mask_1, mask_2 real , dimension(n) :: delta_b !!, us, bs, qs real :: r_crit, sqrt_drag_min @@ -257,7 +288,7 @@ _PURE subroutine monin_obukhov_solve_zeta(error, zeta_min, max_iter, small, & integer, intent(in ) :: n real , intent(in ), dimension(n) :: rich, z, z0, zt, zq logical, intent(in ), dimension(n) :: mask - real , intent( out), dimension(n) :: f_m, f_t, f_q + real , intent( out), dimension(CRAY_DIM) :: f_m, f_t, f_q integer, intent( out) :: ier diff --git a/src/atmos_param/vert_turb_driver/vert_turb_driver.F90 b/src/atmos_param/vert_turb_driver/vert_turb_driver.F90 index 01784c0f3..364cb54d4 100644 --- a/src/atmos_param/vert_turb_driver/vert_turb_driver.F90 +++ b/src/atmos_param/vert_turb_driver/vert_turb_driver.F90 @@ -61,7 +61,7 @@ module vert_turb_driver_mod use fms_mod, only: mpp_pe, mpp_root_pe, stdlog, & error_mesg, open_namelist_file, file_exist, & check_nml_error, close_file, FATAL, & - write_version_number + write_version_number, close_file use field_manager_mod, only: MODEL_ATMOS @@ -589,7 +589,9 @@ subroutine vert_turb_driver_init (lonb, latb, id, jd, kd, axes, Time, & !--------------- read namelist ------------------ #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=vert_turb_driver_nml, iostat=io) + unit = open_namelist_file() + read (unit, nml=vert_turb_driver_nml, iostat=io) + call close_file(unit) ierr = check_nml_error(io,'vert_turb_driver_nml') #else if (file_exist('input.nml')) then diff --git a/src/atmos_shared/vert_advection/vert_advection.F90 b/src/atmos_shared/vert_advection/vert_advection.F90 index 02b34f34e..4a033fec4 100644 --- a/src/atmos_shared/vert_advection/vert_advection.F90 +++ b/src/atmos_shared/vert_advection/vert_advection.F90 @@ -591,7 +591,7 @@ subroutine compute_weights ( dz, zwt ) ! coefficients/weights for computing values at grid box interfaces ! only recompute coefficients for a column when layer depth has changed - + !dir$ IVDEP do j = 1, size(dz,2) do i = 1, size(dz,1) @@ -621,6 +621,7 @@ subroutine compute_weights ( dz, zwt ) enddo dzs(i,j,:) = dz(i,j,:) zwts(0:3,i,j,:) = zwt(0:3,i,j,:) + else ! use previously computed coefficients diff --git a/src/atmos_solo/atmos_model.F90 b/src/atmos_solo/atmos_model.F90 index 56d4e6a56..d43dfa7db 100644 --- a/src/atmos_solo/atmos_model.F90 +++ b/src/atmos_solo/atmos_model.F90 @@ -42,14 +42,16 @@ program atmos_model NOLEAP, NO_CALENDAR, set_calendar_type, & set_date, get_date -use fms_mod, only: file_exist, check_nml_error, & - error_mesg, FATAL, WARNING, & - mpp_pe, mpp_root_pe, fms_init, fms_end, & - stdlog, stdout, write_version_number, & - open_restart_file, & - mpp_clock_id, mpp_clock_begin, & - mpp_clock_end, CLOCK_COMPONENT, set_domain, & - nullify_domain, uppercase +use fms_mod, only: file_exist, check_nml_error, & + error_mesg, FATAL, WARNING, & + mpp_pe, mpp_root_pe, fms_init, fms_end, & + stdlog, stdout, write_version_number, & + open_restart_file, & + mpp_clock_id, mpp_clock_begin, & + mpp_clock_end, CLOCK_COMPONENT, set_domain, & + nullify_domain, uppercase, open_namelist_file,& + close_file + use fms_io_mod, only: fms_io_exit use mpp_mod, only: mpp_set_current_pelist @@ -173,7 +175,9 @@ subroutine atmos_model_init !----- read namelist ------- #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=main_nml, iostat=io) + unit = open_namelist_file () + read (unit, nml=main_nml, iostat=io) + call mpp_close (unit) ierr = check_nml_error(io, 'main_nml') #else unit = open_namelist_file ( ) diff --git a/src/atmos_spectral/driver/solo/atmosphere.F90 b/src/atmos_spectral/driver/solo/atmosphere.F90 index 3091e1fe0..ed7c1b67a 100644 --- a/src/atmos_spectral/driver/solo/atmosphere.F90 +++ b/src/atmos_spectral/driver/solo/atmosphere.F90 @@ -27,8 +27,9 @@ module atmosphere_mod use fms_mod, only: open_namelist_file, close_file #endif -use fms_mod, only: set_domain, write_version_number, field_size, file_exist, stdlog, & - mpp_pe, mpp_root_pe, error_mesg, FATAL, read_data, write_data, nullify_domain +use fms_mod, only: set_domain, write_version_number, field_size, file_exist, stdlog, check_nml_error, & + mpp_pe, mpp_root_pe, error_mesg, FATAL, read_data, write_data, nullify_domain, & + open_namelist_file, close_file use constants_mod, only: grav, pi @@ -121,7 +122,7 @@ subroutine atmosphere_init(Time_init, Time, Time_step_in) type (time_type), intent(in) :: Time_init, Time, Time_step_in -integer :: seconds, days, lon_max, lat_max, ntr, nt, i, j, nml_unit, io, stdlog_unit +integer :: seconds, days, lon_max, lat_max, ntr, nt, i, j, nml_unit, io, stdlog_unit, unit, ierr integer, dimension(4) :: siz real, dimension(2) :: time_pointers character(len=64) :: file, tr_name @@ -132,7 +133,10 @@ subroutine atmosphere_init(Time_init, Time, Time_step_in) call write_version_number(version, tagname) #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=atmosphere_nml, iostat=io) + unit = open_namelist_file ( ) + read (unit, nml=atmosphere_nml, iostat=io) + call close_file(unit) + ierr = check_nml_error(io, 'atmosphere_nml') #else if ( file_exist('input.nml') ) then nml_unit = open_namelist_file() diff --git a/src/atmos_spectral/init/polvani_2007.F90 b/src/atmos_spectral/init/polvani_2007.F90 index 677493dd0..bdff1fead 100644 --- a/src/atmos_spectral/init/polvani_2007.F90 +++ b/src/atmos_spectral/init/polvani_2007.F90 @@ -353,7 +353,7 @@ subroutine compute_LC1 if(k == num_levels+1) then Tr = T0 else - Tr = T0 + lapse/(zt**-alpha + z(k)**-alpha)**(1./alpha) + Tr = T0 + lapse/(zt** (-alpha) + z(k)** (-alpha))**(1./alpha) endif LC1_t(:,k) = Tr + LC1_t(:,k) enddo diff --git a/src/atmos_spectral/init/spectral_initialize_fields.F90 b/src/atmos_spectral/init/spectral_initialize_fields.F90 index fc0d87232..212ad52e0 100644 --- a/src/atmos_spectral/init/spectral_initialize_fields.F90 +++ b/src/atmos_spectral/init/spectral_initialize_fields.F90 @@ -114,7 +114,7 @@ subroutine spectral_initialize_fields(reference_sea_level_press, triang_trunc, i ! initial spectral fields (and spectrally-filtered) grid fields call trans_grid_to_spherical(tg, ts) -call trans_spherical_to_grid(ts, tg) +call trans_spherical_to_grid(ts, tg) call trans_grid_to_spherical(ln_psg, ln_ps) call trans_spherical_to_grid(ln_ps, ln_psg) diff --git a/src/atmos_spectral/init/vert_coordinate.F90 b/src/atmos_spectral/init/vert_coordinate.F90 index 0f6cb10da..2551a0ce4 100644 --- a/src/atmos_spectral/init/vert_coordinate.F90 +++ b/src/atmos_spectral/init/vert_coordinate.F90 @@ -38,7 +38,7 @@ module vert_coordinate_mod use fms_mod, only: mpp_pe, mpp_root_pe, error_mesg, FATAL, & write_version_number, stdlog, close_file, & - check_nml_error + check_nml_error, open_namelist_file, close_file use constants_mod, only: pi @@ -193,7 +193,9 @@ subroutine read_namelist (a, b) bk = 0. #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=vert_coordinate_nml, iostat=io) + namelist_unit = open_namelist_file() + read (namelist_unit, nml=vert_coordinate_nml, iostat=io) + call close_file(namelist_unit) ierr = check_nml_error(io, 'vert_coordinate_nml') #else namelist_unit = open_namelist_file() diff --git a/src/atmos_spectral/model/spectral_dynamics.F90 b/src/atmos_spectral/model/spectral_dynamics.F90 index 0a9b2b310..0820e28d1 100644 --- a/src/atmos_spectral/model/spectral_dynamics.F90 +++ b/src/atmos_spectral/model/spectral_dynamics.F90 @@ -30,7 +30,7 @@ module spectral_dynamics_mod use fms_mod, only: mpp_pe, mpp_root_pe, error_mesg, NOTE, FATAL, write_version_number, stdlog, & close_file, open_restart_file, file_exist, set_domain, & read_data, write_data, check_nml_error, lowercase, uppercase, mpp_npes, & - field_size + field_size, open_namelist_file use constants_mod, only: rdgas, rvgas, grav, cp_air, omega, radius, pi @@ -253,22 +253,23 @@ subroutine spectral_dynamics_init(Time, Time_step_in, tracer_attributes, dry_mod if(module_is_initialized) return #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=spectral_dynamics_nml, iostat=io) - ierr = check_nml_error(io, 'spectral_dynamics_nml') + unit = open_namelist_file ( ) + read (unit, nml=spectral_dynamics_nml, iostat=io) + call close_file(unit) + ierr = check_nml_error(io, 'spectral_dynamics_nml') #else - unit = open_namelist_file() - ierr=1 - do while (ierr /= 0) - read(unit, nml=spectral_dynamics_nml, iostat=io, end=20) - ierr = check_nml_error (io, 'spectral_dynamics_nml') - enddo -20 call close_file (unit) + unit = open_namelist_file() + ierr=1 + do while (ierr /= 0) + read(unit, nml=spectral_dynamics_nml, iostat=io, end=20) + ierr = check_nml_error (io, 'spectral_dynamics_nml') + enddo + 20 call close_file (unit) #endif call write_version_number(version, tagname) if(mpp_pe() == mpp_root_pe()) write (stdlog(), nml=spectral_dynamics_nml) call write_version_number(tracer_type_version, tracer_type_tagname) - Time_step = Time_step_in Alarm_interval = set_time(print_interval(2), print_interval(1)) Alarm_time = Time + Alarm_interval diff --git a/src/atmos_spectral/tools/grid_fourier.F90 b/src/atmos_spectral/tools/grid_fourier.F90 index 821ce5ff4..6af9f4153 100644 --- a/src/atmos_spectral/tools/grid_fourier.F90 +++ b/src/atmos_spectral/tools/grid_fourier.F90 @@ -146,6 +146,7 @@ function trans_grid_to_fourier_3d(grid) result(fourier) grid1(1:size(grid,1),:,:) = grid grid1(size(grid,1)+1,:,:) = grid(1,:,:) + fourier = fft_grid_to_fourier(grid1) return diff --git a/src/atmos_spectral/tools/spherical_fourier.F90 b/src/atmos_spectral/tools/spherical_fourier.F90 index e5396fac5..9f5b0a1a0 100644 --- a/src/atmos_spectral/tools/spherical_fourier.F90 +++ b/src/atmos_spectral/tools/spherical_fourier.F90 @@ -202,7 +202,7 @@ subroutine trans_spherical_to_fourier_3d(spherical,fourier) nj = je-js+1 !this version requires grid domains equal in Y if( debug )write( 0,'(a,i2,2i4)' )'pe, nd, nj=', mpp_pe(), nd, nj - if( mod(ns,2).EQ.0 )then + if( mod(ns,2).EQ.0 ) then neven = ns nodd = ns+1 else @@ -256,7 +256,6 @@ subroutine trans_spherical_to_fourier_3d(spherical,fourier) end if end do end do DOMAIN_LOOP - return end subroutine trans_spherical_to_fourier_3d diff --git a/src/atmos_spectral/tools/spherical_fourier_fftw.F90 b/src/atmos_spectral/tools/spherical_fourier_fftw.F90 new file mode 100644 index 000000000..e7ea97d7a --- /dev/null +++ b/src/atmos_spectral/tools/spherical_fourier_fftw.F90 @@ -0,0 +1,180 @@ +module spherical_fourier_fftw_mod +use fftw3 +use fms_mod, only: mpp_pe, mpp_root_pe, error_mesg, FATAL, write_version_number + +use mpp_mod, only: mpp_error + +use mpp_domains_mod, only: domain1D, mpp_get_compute_domains, mpp_get_domain_components, mpp_get_layout + +use constants_mod, only: pi + +use spec_mpp_mod, only: get_grid_domain, grid_domain, get_spec_domain + +use spherical_mod, only: spherical_init, spherical_end, & + compute_lon_deriv_cos, compute_lat_deriv_cos, & + compute_laplacian, compute_vor, compute_div, & + get_spherical_wave, get_fourier_wave, & + get_eigen_laplacian, compute_gradient_cos, & + compute_ucos_vcos, compute_vor_div, & + triangular_truncation, rhomboidal_truncation + +use gauss_and_legendre_mod, only: compute_legendre, compute_gaussian + +implicit none +private + +character(len=128), parameter :: version = '$Id: spherical_fourier_fftw.F90,v 13.0 2019 fms Exp $' +character(len=128), parameter :: tagname = '$Name: siena_201211 $' + +interface trans_spherical_to_fourier_fftw + module procedure trans_spherical_to_fourier_fftw_3d, & + trans_spherical_to_fourier_fftw_2d +end interface +interface trans_fourier_to_spherical_fftw + module procedure trans_fourier_to_spherical_fftw_3d, & + trans_fourier_to_spherical_fftw_2d +end interface + +interface trans_grid_to_spherical_fftw + module procedure trans_grid_to_spherical_fftw_3d, & + trans_grid_to_spherical_fftw_2d +end interface + +logical :: module_is_initialized = .false. +integer :: is, ie, js, je, ms, me, ns, ne +integer :: fourier_inc, num_fourier, fourier_max, lat_max, & + num_spherical +logical :: make_symmetric_local + +public :: spherical_fourier_fftw_init, trans_grid_to_spherical_fftw + +contains + +!----------------------------------------------------------------------- +subroutine spherical_fourier_fftw_init(radius, lat_max_in, num_fourier_in, & + fourier_inc_in, num_spherical_in, south_to_north, make_symmetric) +!----------------------------------------------------------------------- + +real, intent(in) :: radius +integer, intent(in) :: lat_max_in +integer, intent(in) :: num_fourier_in +integer, intent(in) :: fourier_inc_in +integer, intent(in) :: num_spherical_in +logical, intent(in), optional :: south_to_north, make_symmetric + +call write_version_number(version, tagname) + + +if(present(make_symmetric)) then + make_symmetric_local = make_symmetric +else + make_symmetric_local = .false. +end if + +call get_grid_domain(is, ie, js, je) +call get_spec_domain(ms, me, ns, ne) + +lat_max = lat_max_in + +fourier_inc = fourier_inc_in +num_spherical = num_spherical_in +num_fourier = num_fourier_in +fourier_max = num_fourier*fourier_inc + + + +module_is_initialized = .true. + +return +end subroutine spherical_fourier_fftw_init + + + +!!!! work on this one +subroutine trans_grid_to_spherical_fftw_3d(grid, spherical) + real, intent(inout), dimension (is:,:,:) :: grid + complex, intent(out), dimension (ms:,ns:,:) :: spherical + + ! call forwardsht_real_fortran() + + ! integer :: num_levels + ! type(C_PTR) :: p, plan + ! complex(C_DOUBLE_COMPLEX), pointer :: arr(:,:,:) + + ! num_levels = 25 + + ! p = fftw_alloc_complex(int(ie * je * num_levels, C_SIZE_T)) + ! plan = fftw_plan_dft_r2c_3d(ie, je, num_levels, grid, spherical, 0) + + ! call fftw_execute_dft_r2c(plan, grid, spherical) + +return +end subroutine trans_grid_to_spherical_fftw_3d + +subroutine trans_grid_to_spherical_fftw_2d(grid, spherical) + real, intent(in), dimension (:,:) :: grid + complex, intent (out), dimension (:,:) :: spherical + + + + + + + ! integer :: num_levels + ! type(C_PTR) :: p, plan + ! complex(C_DOUBLE_COMPLEX), pointer :: arr(:,:,:) + + ! num_levels = 25 + + ! p = fftw_alloc_complex(int(ie * je, C_SIZE_T)) + ! plan = fftw_plan_dft_r2c_2d(ie, je, grid, spherical, 0) + + ! call fftw_execute_dft_r2c(plan, grid, spherical) + + + + + +return +end subroutine trans_grid_to_spherical_fftw_2d + + + + + + + + + + + + +subroutine trans_spherical_to_fourier_fftw_3d(spherical, fourier) + complex, intent(in), dimension(ms:,ns:,:) :: spherical + complex, intent(out), dimension(ms:,:,:,:) :: fourier + +end subroutine trans_spherical_to_fourier_fftw_3d + +subroutine trans_spherical_to_fourier_fftw_2d(spherical, fourier) +complex, intent(in), dimension(:,:) :: spherical +complex, intent(out), dimension(:,:,:) :: fourier + +end subroutine trans_spherical_to_fourier_fftw_2d + + + +subroutine trans_fourier_to_spherical_fftw_3d(spherical, fourier) +complex, intent(in), dimension(:,:,:) :: spherical +complex, intent(out), dimension(:,:,:) :: fourier + +end subroutine trans_fourier_to_spherical_fftw_3d + +subroutine trans_fourier_to_spherical_fftw_2d(spherical, fourier) +complex, intent(in), dimension(:,:) :: spherical +complex, intent(out), dimension(:,:,:) :: fourier + +end subroutine trans_fourier_to_spherical_fftw_2d + + + +end module spherical_fourier_fftw_mod \ No newline at end of file diff --git a/src/atmos_spectral/tools/transforms.F90 b/src/atmos_spectral/tools/transforms.F90 index 121580d48..dd9f39e63 100644 --- a/src/atmos_spectral/tools/transforms.F90 +++ b/src/atmos_spectral/tools/transforms.F90 @@ -27,7 +27,7 @@ module transforms_mod use fms_mod, only: open_namelist_file #endif -use fms_mod, only: mpp_pe, mpp_root_pe, error_mesg, FATAL, write_version_number, stdlog, close_file, check_nml_error +use fms_mod, only: mpp_pe, mpp_root_pe, error_mesg, FATAL, write_version_number, stdlog, close_file, check_nml_error, open_namelist_file use mpp_mod, only: mpp_chksum, mpp_error, mpp_npes, mpp_sum, mpp_sync, mpp_sync_self, mpp_transmit @@ -79,7 +79,7 @@ module transforms_mod !--------------------------------------------------------------------------- implicit none -private +! private interface trans_spherical_to_grid module procedure trans_spherical_to_grid_3d, & @@ -246,8 +246,10 @@ subroutine transforms_init(radius, & if(module_is_initialized) return #ifdef INTERNAL_FILE_NML - read (input_nml_file, nml=transforms_nml, iostat=io) - ierr = check_nml_error(io, 'transforms_nml') + namelist_unit = open_namelist_file() + read (namelist_unit, transforms_nml, iostat=io) + call close_file(namelist_unit) + ierr = check_nml_error(io,'transforms_nml') #else namelist_unit = open_namelist_file() ierr=1 @@ -301,7 +303,7 @@ subroutine transforms_init(radius, & ! initialize spherical_fourier (which initializes spherical) call spherical_fourier_init(radius, lat_max, num_fourier, fourier_inc, num_spherical, & - south_to_north=south_to_north_local, make_symmetric=make_symmetric_local) + south_to_north=south_to_north_local, make_symmetric=make_symmetric_local) trunc_fourier = num_fourier @@ -376,12 +378,13 @@ subroutine reset_num_lon_in_transform(num_lon_in, trunc_fourier_in, longitude_or end subroutine reset_num_lon_in_transform !-------------------------------------------------------------------------- - subroutine trans_spherical_to_grid_3d(spherical, grid) + subroutine trans_spherical_to_grid_3d(spherical, grid, block) !-------------------------------------------------------------------------- complex, intent (in), dimension (ms:,ns:,:) :: spherical real, intent(out), dimension (is:,:,:) :: grid real, dimension(num_lon,size(grid,2),size(grid,3)) :: grid_xglobal +logical, intent(in), optional :: block integer(kind=kind(spherical)) :: c1, c2, c3 @@ -390,6 +393,11 @@ subroutine trans_spherical_to_grid_3d(spherical, grid) logical :: grid_x_is_global, spectral_y_is_global type(domain1D) :: spectral_domain_y integer, allocatable :: pelist(:) +logical :: block_comm = .true. + +if(PRESENT(block)) then + block_comm = block +endif if(.not.module_is_initialized) then call error_mesg('trans_spherical_to_grid','transforms module is not initialized', FATAL) @@ -410,6 +418,8 @@ subroutine trans_spherical_to_grid_3d(spherical, grid) if( size(spherical,3).ne.size(grid,3) )& call mpp_error( FATAL, 'TRANS_SPHERICAL_TO_GRID: size(spherical,3).ne.size(grid,3.' ) + + call trans_spherical_to_fourier( spherical, fourier_s ) call mpp_get_compute_domain( spectral_domain, y_is_global=spectral_y_is_global ) if( .NOT.spectral_y_is_global )then @@ -419,7 +429,7 @@ subroutine trans_spherical_to_grid_3d(spherical, grid) call mpp_sum( fourier_s , size(fourier_s(:,:,:,:)), pelist ) end if -call reverse_transpose_fourier( fourier_s, fourier_g ) +call reverse_transpose_fourier( fourier_s, fourier_g, block=block_comm ) fourier_g(trunc_fourier+1:num_lon/2,:,:) = cmplx(0.,0.) call mpp_get_compute_domain( grid_domain, x_is_global=grid_x_is_global ) @@ -442,17 +452,23 @@ subroutine trans_spherical_to_grid_3d(spherical, grid) end subroutine trans_spherical_to_grid_3d !-------------------------------------------------------------------------- - subroutine trans_spherical_to_grid_2d(spherical, grid) + subroutine trans_spherical_to_grid_2d(spherical, grid, block) !-------------------------------------------------------------------------- complex, intent (in), dimension (:,:) :: spherical real, intent(out), dimension (:,:) :: grid +logical, optional :: block real, dimension (size(grid,1), size(grid,2), 1) :: grid_3d complex, dimension(size(spherical,1), size(spherical,2), 1) :: spherical_3d +logical :: block_comm = .true. + +if(PRESENT(block)) then + block_comm = block +endif spherical_3d(:,:,1) = spherical(:,:) -call trans_spherical_to_grid_3d(spherical_3d, grid_3d) +call trans_spherical_to_grid_3d(spherical_3d, grid_3d, block_comm) grid(:,:) = grid_3d(:,:,1) return @@ -472,7 +488,6 @@ subroutine trans_grid_to_spherical_3d(grid, spherical, do_truncation) complex, dimension (ms:me, je-js+1, size(grid,3), grid_layout(2)) :: fourier_s integer(kind=kind(spherical)) :: c1, c2, c3 - if(.not.module_is_initialized) then call error_mesg('trans_grid_to_spherical','transforms module is not initialized', FATAL) end if @@ -967,14 +982,24 @@ subroutine get_grid_boundaries(lon_boundaries, lat_boundaries,global) end subroutine get_grid_boundaries !------------------------------------------------------------------------- -subroutine reverse_transpose_fourier( fourier_s, fourier_g ) +subroutine reverse_transpose_fourier( fourier_s, fourier_g , block) !------------------------------------------------------------------------- complex, intent(in) :: fourier_s(:,:,:,0:) complex, intent(out) :: fourier_g(0:,:,:) + logical, optional, intent(in) :: block complex, dimension(xmaxsize*ymaxsize*size(fourier_s,3)) :: get_data integer :: i,j,k, jj, jp, jm, pp, pm, nput, nget, jpos type(domain1D) :: spectral_domain_x, grid_domain_y integer, dimension(0:grid_layout(2)-1) :: pelist, ygridsize, xspecsize, xsbegin, xsend + real :: start_time = 0.0, stop_time = 0.0 + + logical :: block_comm = .true. + + if(PRESENT(block)) then + block_comm = block + endif + + if(.not.module_is_initialized) then call error_mesg('reverse_transpose_fourier','transforms module is not initialized', FATAL) @@ -987,15 +1012,17 @@ subroutine reverse_transpose_fourier( fourier_s, fourier_g ) call mpp_get_compute_domains( spectral_domain_x, xsbegin, xsend, xspecsize ) nput = size(fourier_s,1)*size(fourier_s,2)*size(fourier_s,3) fourier_g(ms:me,:,:) = fourier_s(:,:,:,jpos) + do jj = 1,grid_layout(2)-1 jp = mod(jpos+jj,grid_layout(2)) jm = mod(jpos-jj+grid_layout(2),grid_layout(2)) pp = pelist(jp) pm = pelist(jm) + nget = xspecsize(jm)*ygridsize(jm)*size(fourier_s,3) ! Force use of "scalar", integer pointer mpp interface call mpp_transmit( put_data=fourier_s(1,1,1,jp), plen=nput, to_pe=pp, & - get_data=get_data(1), glen=nget, from_pe=pm ) + get_data=get_data(1), glen=nget, from_pe=pm,block=block_comm) nget = 0 do k = 1,size(fourier_g,3) do j = 1,size(fourier_g,2) @@ -1006,6 +1033,7 @@ subroutine reverse_transpose_fourier( fourier_s, fourier_g ) end do end do end do + call mpp_sync() return end subroutine reverse_transpose_fourier diff --git a/src/coupler/surface_flux.F90 b/src/coupler/surface_flux.F90 index 5971afebb..7cea8d191 100644 --- a/src/coupler/surface_flux.F90 +++ b/src/coupler/surface_flux.F90 @@ -35,7 +35,7 @@ module surface_flux_mod ! ! ============================================================================ -use fms_mod, only: FATAL, close_file, mpp_pe, mpp_root_pe, write_version_number, error_mesg, FATAL +use fms_mod, only: FATAL, close_file, mpp_pe, mpp_root_pe, write_version_number, error_mesg use fms_mod, only: file_exist, check_nml_error, open_namelist_file, stdlog use monin_obukhov_mod, only: mo_drag, mo_profile use frierson_monin_obukhov_mod, only: frierson_mo_drag diff --git a/src/extra/env/bluepebble b/src/extra/env/bluepebble new file mode 100644 index 000000000..1c365a7e1 --- /dev/null +++ b/src/extra/env/bluepebble @@ -0,0 +1,19 @@ +echo loadmodules for Blue pebble +#module purge + +# Load modules +module load lang/python/anaconda/3.7-2019.03.biopython +module load lang/intel-parallel-studio-xe/2019.u3 +module load lib/netcdf/4.70 +module load tools/git/2.22.0 + +# Compiler variables +export F90=mpiifort +#export F90=mpif90 +export CC=mpiicc + +# Find the NetCDF library +export CPATH="/sw/lib/netcdf-c-4.7.0/include:${CPATH}" + + +export PATH=/usr/bin:$PATH diff --git a/src/extra/env/bristol-bc3 b/src/extra/env/bristol-bc3 new file mode 100644 index 000000000..22d0060d1 --- /dev/null +++ b/src/extra/env/bristol-bc3 @@ -0,0 +1,23 @@ +echo loadmodules for Bristol Bluecrystal 3 + + +# No purge as causes crash + +module add shared default-environment +module load languages/python-anaconda-4.2-3.5 + +module add intel-cluster-studio/compiler/64/13.1/117 +module add openmpi/intel/64/1.6.5 +module add tools/intel_builds/tau-2.23.1-openmpi +module add libraries/intel_builds/netcdf-4.3 +module add tools/git-2.18.0 +module add fftw3/intel/64/3.3.3 + +# source psxevars +export FC=mpif90 +export CC=mpicc + +export CPATH=$CPATH:/cm/shared/apps/fftw/intel/64/3.3.3/include +export CPATH=$CPATH:/cm/shared/libraries/intel_build/netcdf-4/include +export CPATH=$CPATH:/cm/shared/libraries/intel_build/hdf5-1.8.12/include +export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/cm/shared/apps/fftw/intel/64/3.3.3/lib diff --git a/src/extra/env/bristol-bc3-gfortran b/src/extra/env/bristol-bc3-gfortran new file mode 100644 index 000000000..1a0966250 --- /dev/null +++ b/src/extra/env/bristol-bc3-gfortran @@ -0,0 +1,10 @@ +echo Loading basic gfortran environment + +# this defaults to ia64, but we will use gfortran, not ifort +module load openmpi/gcc/64/2.1.1 +module load netcdf/gcc/64/4.2 +module load fftw3/gcc/64/3.3.3 +module load netcdf/gfortran/64/4.2 + +export F90=mpif90 +export CC=mpicc diff --git a/src/extra/env/bristol-bc4 b/src/extra/env/bristol-bc4 index f3c48ada2..0a6492e23 100644 --- a/src/extra/env/bristol-bc4 +++ b/src/extra/env/bristol-bc4 @@ -5,6 +5,8 @@ module purge module load libs/netcdf/4.4.1.1.MPI module load build/gcc-7.2.0 module load tools/git/2.18.0 +module load libs/fftw/3.3.6 + export CPATH=$CPATH:/mnt/storage/software/libraries/intel/netcdf-4.4.1.1-mpi/include diff --git a/src/extra/env/bristol-bc4-gfortran b/src/extra/env/bristol-bc4-gfortran new file mode 100644 index 000000000..be08d63da --- /dev/null +++ b/src/extra/env/bristol-bc4-gfortran @@ -0,0 +1,27 @@ +echo loadmodules for Bristol Bluecrystal 4 + +module purge + +#module load languages/anaconda3 +#module load OpenMPI/3.0.0-GCC-7.2.0-2.29 +#module load libs/netcdf/4.4.1.1.MPI +#module load libs/fftw/3.3.6 + +# GNU compiler +module load languages/anaconda3 +module load gompi/2017a +module load libs/fftw/3.3.6 + +#export CPATH=$CPATH:/mnt/storage/software/libraries/intel/netcdf-4.4.1.1-mpi/include + +export F90=mpifort +export CC=mpicc + + +#echo Loading basic gfortran environment + +# this defaults to ia64, but we will use gfortran, not ifort +#export GFDL_MKMF_TEMPLATE=gfort + +#export F90=mpifort +#export CC=mpicc diff --git a/src/extra/env/cray_aarch64 b/src/extra/env/cray_aarch64 new file mode 100644 index 000000000..c54eec2dd --- /dev/null +++ b/src/extra/env/cray_aarch64 @@ -0,0 +1,6 @@ +echo loadmodules for Isambard +#module purge +module load cray-netcdf-hdf5parallel/4.6.1.3 +module load cray-python/3.6.5.6 +module load PrgEnv-cray/6.0.5 +GFDL_MKMF_TEMPLATE=cray_aarch64 diff --git a/src/extra/env/isambard_armhpc b/src/extra/env/isambard_armhpc new file mode 100644 index 000000000..b098e17e3 --- /dev/null +++ b/src/extra/env/isambard_armhpc @@ -0,0 +1 @@ +echo loadmodules for Isambard_armhpc \ No newline at end of file diff --git a/src/extra/env/isambard_gfortran b/src/extra/env/isambard_gfortran new file mode 100644 index 000000000..6dc39ccff --- /dev/null +++ b/src/extra/env/isambard_gfortran @@ -0,0 +1,18 @@ +echo load modules for Isambard GNU + +# export environment variables +export F90=ftn +export CC=cc +export EXECUTION_TYPE=APRUN +export GFDL_MKMF_TEMPLATE=isambard_gfortran + +# clear environment +#module purge + +# load modules +module load PrgEnv-gnu/6.0.5 +module load gdb4hpc/3.0.10 +module load cray-netcdf-hdf5parallel/4.6.1.3 + +# for performance analysis +module load perftools diff --git a/src/extra/env/isca b/src/extra/env/isca index a4cf1f14c..591dfd61b 100644 --- a/src/extra/env/isca +++ b/src/extra/env/isca @@ -4,7 +4,8 @@ module purge module load intel/2016b module load HDF5/1.8.16-intel-2016b module load netCDF-Fortran/4.4.2-intel-2016b +module load FFTW/3.3.5-intel-2016b export F90=mpiifort export CC=mpiicc -export NETCDF_LIBS = `nf-config --fflags --flibs` \ No newline at end of file +export NETCDF_LIBS = `nf-config --fflags --flibs` diff --git a/src/extra/model/dry_fftw/field_table b/src/extra/model/dry_fftw/field_table new file mode 100644 index 000000000..c2d63b9dd --- /dev/null +++ b/src/extra/model/dry_fftw/field_table @@ -0,0 +1,10 @@ + +"TRACER", "atmos_mod", "sphum" + "longname", "specific humidity" + "units", "kg/kg" + "numerical_representation", "grid" + "hole_filling", "off" + "advect_vert", "finite_volume_parabolic" + "robert_filter", "on" + "profile_type", "fixed", "surface_value=0.0" / + diff --git a/src/extra/model/dry_fftw/path_names b/src/extra/model/dry_fftw/path_names new file mode 100644 index 000000000..287f61a78 --- /dev/null +++ b/src/extra/model/dry_fftw/path_names @@ -0,0 +1,206 @@ +atmos_param/diffusivity/diffusivity.F90 +atmos_param/edt/edt.F90 +atmos_param/entrain/entrain.F90 +atmos_param/hs_forcing/hs_forcing.F90 +atmos_param/lscale_cond/lscale_cond.F90 +atmos_param/my25_turb/my25_turb.F90 +atmos_param/qe_moist_convection/qe_moist_convection.F90 +atmos_param/betts_miller/betts_miller.f90 +atmos_param/sea_esf_rad/null/rad_utilities.F90 +atmos_param/shallow_conv/shallow_conv.F90 +atmos_param/cloud_simple/cloud_simple.F90 +atmos_param/cloud_simple/lcl.F90 +atmos_param/cloud_simple/large_scale_cloud.F90 +atmos_param/cloud_simple/marine_strat_cloud.F90 +atmos_param/cloud_simple/cloud_cover_diags.F90 +atmos_param/cloud_simple/cloud_spookie.F90 +atmos_param/stable_bl_turb/stable_bl_turb.F90 +atmos_param/strat_cloud/null/strat_cloud.F90 +atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +atmos_param/qflux/qflux.f90 +atmos_param/monin_obukhov/monin_obukhov_interfaces.h +atmos_param/monin_obukhov/monin_obukhov_kernel.F90 +atmos_param/monin_obukhov/monin_obukhov.F90 +atmos_param/frierson_monin_obukhov/frierson_monin_obukhov.F90 +atmos_param/dry_convection/dry_convection.f90 +atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 +atmos_param/damping_driver/damping_driver.f90 +atmos_param/mg_drag/mg_drag.f90 +atmos_param/cg_drag/cg_drag.f90 +atmos_param/topo_drag/topo_drag.f90 +atmos_param/vert_diff/vert_diff.F90 +atmos_param/vert_turb_driver/vert_turb_driver.F90 +atmos_param/ras/ras.f90 +atmos_shared/interpolator/interpolator.F90 +atmos_shared/vert_advection/vert_advection.F90 +atmos_solo/atmos_model.F90 +atmos_spectral/driver/solo/atmosphere.F90 +atmos_spectral/driver/solo/idealized_moist_phys.F90 +atmos_spectral/driver/solo/mixed_layer.F90 +atmos_spectral/init/ic_from_external_file.F90 +atmos_spectral/init/jablonowski_2006.F90 +atmos_spectral/init/polvani_2004.F90 +atmos_spectral/init/polvani_2007.F90 +atmos_spectral/init/spectral_init_cond.F90 +atmos_spectral/init/spectral_initialize_fields.F90 +atmos_spectral/init/topog_regularization.F90 +atmos_spectral/init/vert_coordinate.F90 +atmos_spectral/model/every_step_diagnostics.F90 +atmos_spectral/model/fv_advection.F90 +atmos_spectral/model/global_integral.F90 +atmos_spectral/model/implicit.F90 +atmos_spectral/model/leapfrog.F90 +atmos_spectral/model/matrix_invert.F90 +atmos_spectral/model/press_and_geopot.F90 +atmos_spectral/model/spectral_damping.F90 +atmos_spectral/model/spectral_dynamics.F90 +atmos_spectral/model/tracer_type.F90 +atmos_spectral/model/water_borrowing.F90 +atmos_spectral/tools/gauss_and_legendre.F90 +atmos_spectral/tools/grid_fourier.F90 +atmos_spectral/tools/spec_mpp.F90 +atmos_spectral/tools/spherical.F90 +atmos_spectral/tools/spherical_fourier.F90 +atmos_spectral/tools/transforms.F90 +coupler/surface_flux.F90 +shared/axis_utils/axis_utils.F90 +shared/constants/constants.F90 +shared/astronomy/astronomy.f90 +shared/diag_manager/diag_axis.F90 +shared/diag_manager/diag_data.F90 +shared/diag_manager/diag_grid.F90 +shared/diag_manager/diag_manager.F90 +shared/diag_manager/diag_output.F90 +shared/diag_manager/diag_table.F90 +shared/diag_manager/diag_util.F90 +shared/fft/fft99.F90 +shared/fft/fft.F90 +shared/field_manager/field_manager.F90 +shared/field_manager/fm_util.F90 +shared/field_manager/parse.inc +shared/fms/fms.F90 +shared/fms/fms_io.F90 +shared/fms/read_data_2d.inc +shared/fms/read_data_3d.inc +shared/fms/read_data_4d.inc +shared/fms/test_fms_io.F90 +shared/fms/write_data.inc +shared/horiz_interp/horiz_interp_bicubic.F90 +shared/horiz_interp/horiz_interp_bilinear.F90 +shared/horiz_interp/horiz_interp_conserve.F90 +shared/horiz_interp/horiz_interp.F90 +shared/horiz_interp/horiz_interp_spherical.F90 +shared/horiz_interp/horiz_interp_type.F90 +shared/include/fms_platform.h +shared/memutils/memuse.c +shared/memutils/memutils.F90 +shared/mosaic/constant.h +shared/mosaic/create_xgrid.c +shared/mosaic/create_xgrid.h +shared/mosaic/gradient_c2l.c +shared/mosaic/gradient_c2l.h +shared/mosaic/gradient.F90 +shared/mosaic/grid.F90 +shared/mosaic/interp.c +shared/mosaic/interp.h +shared/mosaic/mosaic.F90 +shared/mosaic/mosaic_util.c +shared/mosaic/mosaic_util.h +shared/mosaic/read_mosaic.c +shared/mosaic/read_mosaic.h +shared/mpp/affinity.c +shared/mpp/include/mpp_chksum.h +shared/mpp/include/mpp_chksum_int.h +shared/mpp/include/mpp_chksum_scalar.h +shared/mpp/include/mpp_comm.inc +shared/mpp/include/mpp_comm_mpi.inc +shared/mpp/include/mpp_comm_nocomm.inc +shared/mpp/include/mpp_comm_sma.inc +shared/mpp/include/mpp_data_mpi.inc +shared/mpp/include/mpp_data_nocomm.inc +shared/mpp/include/mpp_data_sma.inc +shared/mpp/include/mpp_define_nest_domains.inc +shared/mpp/include/mpp_do_check.h +shared/mpp/include/mpp_do_checkV.h +shared/mpp/include/mpp_do_get_boundary.h +shared/mpp/include/mpp_do_global_field.h +shared/mpp/include/mpp_domains_comm.inc +shared/mpp/include/mpp_domains_define.inc +shared/mpp/include/mpp_domains_misc.inc +shared/mpp/include/mpp_domains_reduce.inc +shared/mpp/include/mpp_domains_util.inc +shared/mpp/include/mpp_do_redistribute.h +shared/mpp/include/mpp_do_update_ad.h +shared/mpp/include/mpp_do_update.h +shared/mpp/include/mpp_do_update_nest.h +shared/mpp/include/mpp_do_update_nonblock.h +shared/mpp/include/mpp_do_updateV_ad.h +shared/mpp/include/mpp_do_updateV.h +shared/mpp/include/mpp_do_updateV_nonblock.h +shared/mpp/include/mpp_error_a_a.h +shared/mpp/include/mpp_error_a_s.h +shared/mpp/include/mpp_error_s_a.h +shared/mpp/include/mpp_error_s_s.h +shared/mpp/include/mpp_gather.h +shared/mpp/include/mpp_get_boundary.h +shared/mpp/include/mpp_global_field.h +shared/mpp/include/mpp_global_reduce.h +shared/mpp/include/mpp_global_sum_ad.h +shared/mpp/include/mpp_global_sum.h +shared/mpp/include/mpp_global_sum_tl.h +shared/mpp/include/mpp_io_connect.inc +shared/mpp/include/mpp_io_misc.inc +shared/mpp/include/mpp_io_read.inc +shared/mpp/include/mpp_io_util.inc +shared/mpp/include/mpp_io_write.inc +shared/mpp/include/mpp_read_2Ddecomp.h +shared/mpp/include/mpp_reduce_mpi.h +shared/mpp/include/mpp_reduce_nocomm.h +shared/mpp/include/mpp_reduce_sma.h +shared/mpp/include/mpp_sum.inc +shared/mpp/include/mpp_sum_mpi.h +shared/mpp/include/mpp_sum_nocomm.h +shared/mpp/include/mpp_sum_sma.h +shared/mpp/include/mpp_transmit.inc +shared/mpp/include/mpp_transmit_mpi.h +shared/mpp/include/mpp_transmit_nocomm.h +shared/mpp/include/mpp_transmit_sma.h +shared/mpp/include/mpp_update_domains2D_ad.h +shared/mpp/include/mpp_update_domains2D.h +shared/mpp/include/mpp_update_domains2D_nonblock.h +shared/mpp/include/mpp_update_nest_domains.h +shared/mpp/include/mpp_util.inc +shared/mpp/include/mpp_util_mpi.inc +shared/mpp/include/mpp_util_nocomm.inc +shared/mpp/include/mpp_util_sma.inc +shared/mpp/include/mpp_write_2Ddecomp.h +shared/mpp/include/mpp_write.h +shared/mpp/include/system_clock.h +shared/mpp/mpp_data.F90 +shared/mpp/mpp_domains.F90 +shared/mpp/mpp.F90 +shared/mpp/mpp_io.F90 +shared/mpp/mpp_memutils.F90 +shared/mpp/mpp_parameter.F90 +shared/mpp/mpp_pset.F90 +shared/mpp/mpp_utilities.F90 +shared/mpp/nsclock.c +shared/mpp/test_mpp_domains.F90 +shared/mpp/test_mpp.F90 +shared/mpp/test_mpp_io.F90 +shared/mpp/test_mpp_pset.F90 +shared/mpp/threadloc.c +shared/platform/platform.F90 +shared/random_numbers/MersenneTwister.F90 +shared/random_numbers/random_numbers.F90 +shared/sat_vapor_pres/sat_vapor_pres.F90 +shared/sat_vapor_pres/sat_vapor_pres_k.F90 +shared/time_interp/time_interp_external.F90 +shared/time_interp/time_interp.F90 +shared/time_manager/get_cal_time.F90 +shared/time_manager/time_manager.F90 +shared/topography/gaussian_topog.F90 +shared/topography/topography.F90 +shared/tracer_manager/tracer_manager.F90 +shared/tridiagonal/tridiagonal.F90 +shared/fft/fftw.F90 diff --git a/src/extra/model/grey_fftw/field_table b/src/extra/model/grey_fftw/field_table new file mode 100644 index 000000000..c2d63b9dd --- /dev/null +++ b/src/extra/model/grey_fftw/field_table @@ -0,0 +1,10 @@ + +"TRACER", "atmos_mod", "sphum" + "longname", "specific humidity" + "units", "kg/kg" + "numerical_representation", "grid" + "hole_filling", "off" + "advect_vert", "finite_volume_parabolic" + "robert_filter", "on" + "profile_type", "fixed", "surface_value=0.0" / + diff --git a/src/extra/model/grey_fftw/path_names b/src/extra/model/grey_fftw/path_names new file mode 100644 index 000000000..094f5ea4d --- /dev/null +++ b/src/extra/model/grey_fftw/path_names @@ -0,0 +1,206 @@ +atmos_param/diffusivity/diffusivity.F90 +atmos_param/edt/edt.F90 +atmos_param/entrain/entrain.F90 +atmos_param/hs_forcing/hs_forcing.F90 +atmos_param/lscale_cond/lscale_cond.F90 +atmos_param/my25_turb/my25_turb.F90 +atmos_param/qe_moist_convection/qe_moist_convection.F90 +atmos_param/betts_miller/betts_miller.f90 +atmos_param/sea_esf_rad/null/rad_utilities.F90 +atmos_param/shallow_conv/shallow_conv.F90 +atmos_param/stable_bl_turb/stable_bl_turb.F90 +atmos_param/strat_cloud/null/strat_cloud.F90 +atmos_param/cloud_simple/cloud_simple.F90 +atmos_param/cloud_simple/lcl.F90 +atmos_param/cloud_simple/large_scale_cloud.F90 +atmos_param/cloud_simple/marine_strat_cloud.F90 +atmos_param/cloud_simple/cloud_cover_diags.F90 +atmos_param/cloud_simple/cloud_spookie.F90 +atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +atmos_param/qflux/qflux.f90 +atmos_param/monin_obukhov/monin_obukhov_interfaces.h +atmos_param/monin_obukhov/monin_obukhov_kernel.F90 +atmos_param/monin_obukhov/monin_obukhov.F90 +atmos_param/frierson_monin_obukhov/frierson_monin_obukhov.F90 +atmos_param/dry_convection/dry_convection.f90 +atmos_param/ras/ras.f90 +atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 +atmos_param/damping_driver/damping_driver.f90 +atmos_param/mg_drag/mg_drag.f90 +atmos_param/cg_drag/cg_drag.f90 +atmos_param/topo_drag/topo_drag.f90 +atmos_param/vert_diff/vert_diff.F90 +atmos_param/vert_turb_driver/vert_turb_driver.F90 +atmos_shared/interpolator/interpolator.F90 +atmos_shared/vert_advection/vert_advection.F90 +atmos_solo/atmos_model.F90 +atmos_spectral/driver/solo/atmosphere.F90 +atmos_spectral/driver/solo/idealized_moist_phys.F90 +atmos_spectral/driver/solo/mixed_layer.F90 +atmos_spectral/init/ic_from_external_file.F90 +atmos_spectral/init/jablonowski_2006.F90 +atmos_spectral/init/polvani_2004.F90 +atmos_spectral/init/polvani_2007.F90 +atmos_spectral/init/spectral_init_cond.F90 +atmos_spectral/init/spectral_initialize_fields.F90 +atmos_spectral/init/topog_regularization.F90 +atmos_spectral/init/vert_coordinate.F90 +atmos_spectral/model/every_step_diagnostics.F90 +atmos_spectral/model/fv_advection.F90 +atmos_spectral/model/global_integral.F90 +atmos_spectral/model/implicit.F90 +atmos_spectral/model/leapfrog.F90 +atmos_spectral/model/matrix_invert.F90 +atmos_spectral/model/press_and_geopot.F90 +atmos_spectral/model/spectral_damping.F90 +atmos_spectral/model/spectral_dynamics.F90 +atmos_spectral/model/tracer_type.F90 +atmos_spectral/model/water_borrowing.F90 +atmos_spectral/tools/gauss_and_legendre.F90 +atmos_spectral/tools/grid_fourier.F90 +atmos_spectral/tools/spec_mpp.F90 +atmos_spectral/tools/spherical.F90 +atmos_spectral/tools/spherical_fourier.F90 +atmos_spectral/tools/transforms.F90 +coupler/surface_flux.F90 +shared/axis_utils/axis_utils.F90 +shared/constants/constants.F90 +shared/astronomy/astronomy.f90 +shared/diag_manager/diag_axis.F90 +shared/diag_manager/diag_data.F90 +shared/diag_manager/diag_grid.F90 +shared/diag_manager/diag_manager.F90 +shared/diag_manager/diag_output.F90 +shared/diag_manager/diag_table.F90 +shared/diag_manager/diag_util.F90 +shared/fft/fft99.F90 +shared/fft/fft.F90 +shared/field_manager/field_manager.F90 +shared/field_manager/fm_util.F90 +shared/field_manager/parse.inc +shared/fms/fms.F90 +shared/fms/fms_io.F90 +shared/fms/read_data_2d.inc +shared/fms/read_data_3d.inc +shared/fms/read_data_4d.inc +shared/fms/test_fms_io.F90 +shared/fms/write_data.inc +shared/horiz_interp/horiz_interp_bicubic.F90 +shared/horiz_interp/horiz_interp_bilinear.F90 +shared/horiz_interp/horiz_interp_conserve.F90 +shared/horiz_interp/horiz_interp.F90 +shared/horiz_interp/horiz_interp_spherical.F90 +shared/horiz_interp/horiz_interp_type.F90 +shared/include/fms_platform.h +shared/memutils/memuse.c +shared/memutils/memutils.F90 +shared/mosaic/constant.h +shared/mosaic/create_xgrid.c +shared/mosaic/create_xgrid.h +shared/mosaic/gradient_c2l.c +shared/mosaic/gradient_c2l.h +shared/mosaic/gradient.F90 +shared/mosaic/grid.F90 +shared/mosaic/interp.c +shared/mosaic/interp.h +shared/mosaic/mosaic.F90 +shared/mosaic/mosaic_util.c +shared/mosaic/mosaic_util.h +shared/mosaic/read_mosaic.c +shared/mosaic/read_mosaic.h +shared/mpp/affinity.c +shared/mpp/include/mpp_chksum.h +shared/mpp/include/mpp_chksum_int.h +shared/mpp/include/mpp_chksum_scalar.h +shared/mpp/include/mpp_comm.inc +shared/mpp/include/mpp_comm_mpi.inc +shared/mpp/include/mpp_comm_nocomm.inc +shared/mpp/include/mpp_comm_sma.inc +shared/mpp/include/mpp_data_mpi.inc +shared/mpp/include/mpp_data_nocomm.inc +shared/mpp/include/mpp_data_sma.inc +shared/mpp/include/mpp_define_nest_domains.inc +shared/mpp/include/mpp_do_check.h +shared/mpp/include/mpp_do_checkV.h +shared/mpp/include/mpp_do_get_boundary.h +shared/mpp/include/mpp_do_global_field.h +shared/mpp/include/mpp_domains_comm.inc +shared/mpp/include/mpp_domains_define.inc +shared/mpp/include/mpp_domains_misc.inc +shared/mpp/include/mpp_domains_reduce.inc +shared/mpp/include/mpp_domains_util.inc +shared/mpp/include/mpp_do_redistribute.h +shared/mpp/include/mpp_do_update_ad.h +shared/mpp/include/mpp_do_update.h +shared/mpp/include/mpp_do_update_nest.h +shared/mpp/include/mpp_do_update_nonblock.h +shared/mpp/include/mpp_do_updateV_ad.h +shared/mpp/include/mpp_do_updateV.h +shared/mpp/include/mpp_do_updateV_nonblock.h +shared/mpp/include/mpp_error_a_a.h +shared/mpp/include/mpp_error_a_s.h +shared/mpp/include/mpp_error_s_a.h +shared/mpp/include/mpp_error_s_s.h +shared/mpp/include/mpp_gather.h +shared/mpp/include/mpp_get_boundary.h +shared/mpp/include/mpp_global_field.h +shared/mpp/include/mpp_global_reduce.h +shared/mpp/include/mpp_global_sum_ad.h +shared/mpp/include/mpp_global_sum.h +shared/mpp/include/mpp_global_sum_tl.h +shared/mpp/include/mpp_io_connect.inc +shared/mpp/include/mpp_io_misc.inc +shared/mpp/include/mpp_io_read.inc +shared/mpp/include/mpp_io_util.inc +shared/mpp/include/mpp_io_write.inc +shared/mpp/include/mpp_read_2Ddecomp.h +shared/mpp/include/mpp_reduce_mpi.h +shared/mpp/include/mpp_reduce_nocomm.h +shared/mpp/include/mpp_reduce_sma.h +shared/mpp/include/mpp_sum.inc +shared/mpp/include/mpp_sum_mpi.h +shared/mpp/include/mpp_sum_nocomm.h +shared/mpp/include/mpp_sum_sma.h +shared/mpp/include/mpp_transmit.inc +shared/mpp/include/mpp_transmit_mpi.h +shared/mpp/include/mpp_transmit_nocomm.h +shared/mpp/include/mpp_transmit_sma.h +shared/mpp/include/mpp_update_domains2D_ad.h +shared/mpp/include/mpp_update_domains2D.h +shared/mpp/include/mpp_update_domains2D_nonblock.h +shared/mpp/include/mpp_update_nest_domains.h +shared/mpp/include/mpp_util.inc +shared/mpp/include/mpp_util_mpi.inc +shared/mpp/include/mpp_util_nocomm.inc +shared/mpp/include/mpp_util_sma.inc +shared/mpp/include/mpp_write_2Ddecomp.h +shared/mpp/include/mpp_write.h +shared/mpp/include/system_clock.h +shared/mpp/mpp_data.F90 +shared/mpp/mpp_domains.F90 +shared/mpp/mpp.F90 +shared/mpp/mpp_io.F90 +shared/mpp/mpp_memutils.F90 +shared/mpp/mpp_parameter.F90 +shared/mpp/mpp_pset.F90 +shared/mpp/mpp_utilities.F90 +shared/mpp/nsclock.c +shared/mpp/test_mpp_domains.F90 +shared/mpp/test_mpp.F90 +shared/mpp/test_mpp_io.F90 +shared/mpp/test_mpp_pset.F90 +shared/mpp/threadloc.c +shared/platform/platform.F90 +shared/random_numbers/MersenneTwister.F90 +shared/random_numbers/random_numbers.F90 +shared/sat_vapor_pres/sat_vapor_pres.F90 +shared/sat_vapor_pres/sat_vapor_pres_k.F90 +shared/time_interp/time_interp_external.F90 +shared/time_interp/time_interp.F90 +shared/time_manager/get_cal_time.F90 +shared/time_manager/time_manager.F90 +shared/topography/gaussian_topog.F90 +shared/topography/topography.F90 +shared/tracer_manager/tracer_manager.F90 +shared/tridiagonal/tridiagonal.F90 +shared/fft/fftw.F90 diff --git a/src/extra/python/isca/__init__.py b/src/extra/python/isca/__init__.py index 830b10774..1cdf81c64 100644 --- a/src/extra/python/isca/__init__.py +++ b/src/extra/python/isca/__init__.py @@ -83,4 +83,4 @@ def emit(self, event, *args, **kwargs): from isca.experiment import Experiment, DiagTable, Namelist, FailedRunError -from isca.codebase import IscaCodeBase, SocratesCodeBase, DryCodeBase, GreyCodeBase, ShallowCodeBase, BarotropicCodeBase, ColumnCodeBase, SocColumnCodeBase +from isca.codebase import IscaCodeBase, SocratesCodeBase, DryCodeBase, GreyCodeBase, ShallowCodeBase, BarotropicCodeBase, ColumnCodeBase, SocColumnCodeBase, DryCodeBaseFFTW diff --git a/src/extra/python/isca/codebase.py b/src/extra/python/isca/codebase.py index c411ff928..f79cf906c 100644 --- a/src/extra/python/isca/codebase.py +++ b/src/extra/python/isca/codebase.py @@ -119,7 +119,9 @@ def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, # read path names from the default file self.path_names = [] + self.path_names_to_append = [] #Adding option so that individual path names can be appended to existing lists self.compile_flags = [] # users can append to this to add additional compiler options + self.precision_compile_flags = ['-DOVERLOAD_C8'] #Default is to use double precision. User can change to single precision using cb.use_single_precision() before compile step @property def code_is_available(self): @@ -255,11 +257,13 @@ def compile(self, debug=False, optimisation=None): # compile_flags.append('-O%d' % optimisation) compile_flags.extend(self.compile_flags) + compile_flags.extend(self.precision_compile_flags) compile_flags_str = ' '.join(compile_flags) # get path_names from the directory if not self.path_names: self.path_names = self.read_path_names(P(self.srcdir, 'extra', 'model', self.name, 'path_names')) + self.path_names.extend(self.path_names_to_append) self.write_path_names(self.path_names) path_names_str = P(self.builddir, 'path_names') @@ -282,6 +286,20 @@ def compile(self, debug=False, optimisation=None): self.log.info('Compilation complete.') + def enable_fftw3(self): + self.log.info('Going to use the FFTW3 options to run the spectral core') + self.compile_flags.append('-DFFTW3') + self.path_names_to_append.append('shared/fft/fftw.F90') + self.executable_name = self.executable_name.strip('.x')+'_fftw.x' + self.builddir = P(self.workdir, 'build', self.executable_name.split('.')[0]) + + def use_single_precision(self): + self.log.info('Going to use single_precision') + self.precision_compile_flags = ['-DOVERLOAD_C4', '-DOVERLOAD_R4'] + self.executable_name = self.executable_name.strip('.x')+'_single.x' + self.builddir = P(self.workdir, 'build', self.executable_name.split('.')[0]) + + class IscaCodeBase(CodeBase): @@ -341,12 +359,18 @@ def simlink_to_soc_code(self): elif GFDL_SOC is None: error_mesg = 'Socrates code is required for SocratesCodebase, but source code is not provided in location GFDL_SOC='+ str(GFDL_SOC) self.log.error(error_mesg) - raise OSError(error_mesg) + raise OSError(error_mesg) + + def run_socrates_single_precision(self): + self.compile_flags.append('-DSOC_SINGLE_PRECISION') + self.log.info('Socrates running with single precision') + self.executable_name = self.executable_name.strip('.x')+'_soc_single.x' + self.builddir = P(self.workdir, 'build', self.executable_name.split('.')[0]) def __init__(self, *args, **kwargs): super(SocratesCodeBase, self).__init__(*args, **kwargs) self.disable_rrtm() - self.simlink_to_soc_code() + self.simlink_to_soc_code() class SocColumnCodeBase(CodeBase): """Isca without RRTM but with the Met Office radiation scheme, Socrates. THIS VERSION FOR SINGLE COLUMN USE. @@ -462,7 +486,11 @@ class DryCodeBase(GreyCodeBase): name = 'dry' executable_name = 'held_suarez.x' - +class DryCodeBaseFFTW(GreyCodeBase): + """The Held-Suarez model using FFTW + """ + name = 'dry_fftw' + executable_name = 'held_suarez.x' class ShallowCodeBase(CodeBase): """The Shallow Water Equations. diff --git a/src/extra/python/isca/experiment.py b/src/extra/python/isca/experiment.py index 9495bd178..b28df5538 100755 --- a/src/extra/python/isca/experiment.py +++ b/src/extra/python/isca/experiment.py @@ -27,6 +27,13 @@ class Experiment(Logger, EventEmitter): """A basic GFDL experiment""" RESOLUTIONS = { + 'T213': { + 'lon_max': 1024, + 'lat_max': 320, + 'num_fourier': 213, + 'num_spherical': 214 + }, + 'T170': { 'lon_max': 512, 'lat_max': 256, @@ -232,8 +239,15 @@ def run(self, i, restart_file=None, use_restart=True, multi_node=False, num_core sh.cp([filename, P(indir, os.path.split(filename)[1])]) if multi_node: - mpirun_opts += ' -bootstrap pbsdsh -f $PBS_NODEFILE' - + try: + nodefile = os.environ['PBS_NODEFILE'] + mpirun_opts += f' -bootstrap pbsdsh -f {nodefile}' + except: + nodefile = os.environ['GFDL_JOB_NODEFILE'] + mpirun_opts += f' -bootstrap slurm -f {nodefile}' + + self.log.info(f'using nodefile {nodefile}, mpirun_opts = {mpirun_opts}') + if use_restart and not restart_file and i == 1: # no restart file specified, but we are at first run number self.log.warn('use_restart=True, but restart_file not specified. As this is run 1, assuming spin-up from namelist stated initial conditions so continuing.') @@ -397,4 +411,3 @@ def derive(self, new_experiment_name): # class RunSpec(Logger): # def __init__(self, exp): # self.exp = exp - diff --git a/src/extra/python/isca/templates/compile.sh b/src/extra/python/isca/templates/compile.sh index db35c1384..18d64c878 100755 --- a/src/extra/python/isca/templates/compile.sh +++ b/src/extra/python/isca/templates/compile.sh @@ -63,13 +63,13 @@ if [ $debug == True ]; then echo "Compiling in debug mode" # execute mkmf to create makefile -cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 {{compile_flags}}" +cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML {{compile_flags}}" $mkmf -a $sourcedir -t $template_debug -p $executable -c "$cppDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include else # execute mkmf to create makefile -cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 ${CDEFS} {{compile_flags}}" +cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML ${CDEFS} {{compile_flags}}" $mkmf -a $sourcedir -t $template -p $executable -c "$cppDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include fi diff --git a/src/extra/python/isca/templates/mkmf.template.ia64 b/src/extra/python/isca/templates/mkmf.template.ia64 index dcb29e579..eac0288b4 100644 --- a/src/extra/python/isca/templates/mkmf.template.ia64 +++ b/src/extra/python/isca/templates/mkmf.template.ia64 @@ -18,11 +18,14 @@ NETCDF_LIBS = `nc-config --libs` # -diag-disable 6843: # This suppresses the warning: `warning #6843: A dummy argument with an explicit INTENT(OUT) declaration is not given an explicit value.` of which # there are a lot of instances in the GFDL codebase. + FFLAGS = $(CPPFLAGS) -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-intel -i4 -r8 -g -O2 -diag-disable 6843 -mcmodel large #FFLAGS = $(CPPFLAGS) -fltconsistency -stack_temps -safe_cray_ptr -ftz -shared-intel -assume byterecl -g -O0 -i4 -r8 -check -warn -warn noerrors -debug variable_locations -inline_debug_info -traceback FC = $(F90) LD = $(F90) $(NETCDF_LIBS) #CC = mpicc -LDFLAGS = -lnetcdff -lnetcdf -lmpi -shared-intel + +LDFLAGS = -lnetcdff -lnetcdf -lmpi -shared-intel -lfftw3 -lfftw3f -L/cm/shared/apps/fftw/intel/64/3.3.3/lib + CFLAGS = -D__IFC diff --git a/src/extra/python/isca/templates/mkmf.template.ia64_trace b/src/extra/python/isca/templates/mkmf.template.ia64_trace new file mode 100644 index 000000000..eb6c612e2 --- /dev/null +++ b/src/extra/python/isca/templates/mkmf.template.ia64_trace @@ -0,0 +1,28 @@ +# template for the Intel fortran compiler +# typical use with mkmf +# mkmf -t template.ifc -c"-Duse_libMPI -Duse_netCDF" path_names /usr/local/include +CPPFLAGS = -I/usr/local/include +NETCDF_LIBS = `nc-config --libs` + +# FFLAGS: +# -fpp: Use the fortran preprocessor +# -stack_temps: Put temporary runtime arrays on the stack, not heap. +# -safe_cray_ptr: Cray pointers don't alias other variables. +# -ftz: Denormal numbers are flushed to zero. +# -assume byterecl: Specifies the units for the OPEN statement as bytes. +# -shared-intel: Load intel libraries dynamically +# -i4: 4 byte integers +# -r8: 8 byte reals +# -g: Generate symbolic debugging info in code +# -O2: Level 2 speed optimisations +# -diag-disable 6843: +# This suppresses the warning: `warning #6843: A dummy argument with an explicit INTENT(OUT) declaration is not given an explicit value.` of which +# there are a lot of instances in the GFDL codebase. +FFLAGS = $(CPPFLAGS) -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-intel -i4 -r8 -g -O2 -diag-disable 6843 -trace +#FFLAGS = $(CPPFLAGS) -fltconsistency -stack_temps -safe_cray_ptr -ftz -shared-intel -assume byterecl -g -O0 -i4 -r8 -check -warn -warn noerrors -debug variable_locations -inline_debug_info -traceback +FC = $(F90) +LD = $(F90) $(NETCDF_LIBS) +#CC = mpicc + +LDFLAGS = -lnetcdff -lnetcdf -lmpi -shared-intel +CFLAGS = -D__IFC diff --git a/src/extra/python/isca/templates/mkmf.template.isambard_armhpc b/src/extra/python/isca/templates/mkmf.template.isambard_armhpc new file mode 100644 index 000000000..699691cd5 --- /dev/null +++ b/src/extra/python/isca/templates/mkmf.template.isambard_armhpc @@ -0,0 +1,27 @@ +# template for the Cray Fortran compiler on a Cavium ThunderX2 system +# typical use with mkmf +# mkmf -t template.ifc -c"-Duse_libMPI -Duse_netCDF" path_names /usr/local/include +CPPFLAGS = -I/usr/local/include +NETCDF_LIBS = `nc-config --libs` + +# FFLAGS: +# -g: Generate symbolic debugging info in code +# -s real64: Use 64 bit (8 byte) real numbers +# -ef: Generate module files in lowercase +# -eZ: Use a Fortran preprocessor for macros +# -e0: Initialise all numbers as 0 +# -h develop: Developer debugging for faster compilation +# -dynamic: Use dynamic linking +# -h cpu=arm-thunderx2: Compile for the thunderx2 processor + +#FFLAGS = $(CPPFLAGS) -g -cpp -r8 -O2 -fstack-arrays -mcpu=native +FFLAGS = $(CPPFLAGS) -g -cpp -r8 -O2 -mcpu=native -ffp-contract=fast -ffree-form -ffree-line-length-none -D_NOGETPID +#FFLAGS = $(CPPFLAGS) -fltconsistency -stack_temps -safe_cray_ptr -ftz -shared-intel -assume byterecl -g -O0 -i4 -r8 -check -warn -warn noerrors -debug variable_locations -inline_debug_info -traceback +FC = $(F90) +LD = $(F90) $(NETCDF_LIBS) +#CC = mpicc + +LDFLAGS = -lnetcdff -lnetcdf -g +CFLAGS = -D__IFC + +MONIN_OBUKHOV_KERNEL_FLAGS = -ev \ No newline at end of file diff --git a/src/extra/python/isca/templates/mkmf.template.isambard_cce b/src/extra/python/isca/templates/mkmf.template.isambard_cce new file mode 100644 index 000000000..1ed9ff9cc --- /dev/null +++ b/src/extra/python/isca/templates/mkmf.template.isambard_cce @@ -0,0 +1,24 @@ +# template for the Cray Fortran compiler on a Cavium ThunderX2 system +# typical use with mkmf +# mkmf -t template.ifc -c"-Duse_libMPI -Duse_netCDF" path_names /usr/local/include +CPPFLAGS = -I/usr/local/include +NETCDF_LIBS = `nc-config --libs` + +# FFLAGS: +# -g: Generate symbolic debugging info in code +# -s real64: Use 64 bit (8 byte) real numbers +# -ef: Generate module files in lowercase +# -eZ: Use a Fortran preprocessor for macros +# -e0: Initialise all numbers as 0 +# -dynamic: Use dynamic linking +# -h cpu=arm-thunderx2: Compile for the thunderx2 processor + +FFLAGS = $(CPPFLAGS) -g -s real64 -efZ0 -h develop -dynamic -hsystem_alloc -h cpu=arm-thunderx2 + +FC = $(F90) +LD = $(F90) $(NETCDF_LIBS) + +LDFLAGS = -lnetcdff -lnetcdf -g +CFLAGS = -D__IFC + +MONIN_OBUKHOV_KERNEL_FLAGS = -ev \ No newline at end of file diff --git a/src/extra/python/isca/templates/mkmf.template.isambard_gfortran b/src/extra/python/isca/templates/mkmf.template.isambard_gfortran new file mode 100644 index 000000000..63e38ceb9 --- /dev/null +++ b/src/extra/python/isca/templates/mkmf.template.isambard_gfortran @@ -0,0 +1,22 @@ +# template for the gfortran compiler +# typical use with mkmf +# mkmf -t template.ifc -c"-Duse_libMPI -Duse_netCDF" path_names /usr/local/include +CPPFLAGS = -I/usr/local/include +NETCDF_LIBS = `nf-config --fflags --flibs` + +# FFLAGS: +# -cpp: Use the fortran preprocessor +# -ffree-line-length-none -fno-range-check: Allow arbitrarily long lines +# -fcray-pointer: Cray pointers don't alias other variables +# -fdefault-real-8: 8 byte reals (compatability for some parts of GFDL code) +# -fdefault-double-8: 8 byte doubles (compat. with RRTM) +# -ftree-vectorize: Enable SIMD instructions +# -O2: Level 2 speed optimisations + +FFLAGS = $(CPPFLAGS) $(NETCDF_LIBS) -cpp -fcray-pointer -O2 -ffree-line-length-none -fno-range-check -target-cpu=arm-thunderx2 -ftree-vectorize -march=armv8.1-a -mcpu=thunderx2t99 -mtune=thunderx2t99 -fdefault-real-8 -fdefault-double-8 + +FC = $(F90) +LD = $(F90) $(NETCDF_LIBS) + +LDFLAGS = -lnetcdff -lnetcdf +CFLAGS = -D__IFC diff --git a/src/extra/python/isca/templates/run.sh b/src/extra/python/isca/templates/run.sh index e4a8d606c..77c2c9774 100644 --- a/src/extra/python/isca/templates/run.sh +++ b/src/extra/python/isca/templates/run.sh @@ -19,7 +19,12 @@ if [ $debug == True ]; then echo "Opening idb for debugging" exec idb -gdb {{ executable}} else - exec nice -{{nice_score}} mpirun {{mpirun_opts}} -np {{ num_cores }} {{ execdir }}/{{ executable }} + # defaults to mpirun -> export EXECUTION_TYPE=APRUN for Cray systems + if [[ -z "${EXECUTION_TYPE}" ]] || [ "${EXECUTION_TYPE^^}" = "MPIRUN" ]; then + exec nice -{{nice_score}} mpirun {{mpirun_opts}} -n {{ num_cores }} {{ execdir }}/{{ executable }} + elif [ "${EXECUTION_TYPE^^}" = "APRUN" ]; then + exec nice -{{nice_score}} aprun {{mpirun_opts}} -n {{ num_cores }} {{ execdir }}/{{ executable }} + fi fi err_code=$? diff --git a/src/extra/python/scripts/modified_time_script.py b/src/extra/python/scripts/modified_time_script.py index d4b5c6249..83af3e3e2 100644 --- a/src/extra/python/scripts/modified_time_script.py +++ b/src/extra/python/scripts/modified_time_script.py @@ -6,7 +6,7 @@ from datetime import datetime -def calculate_month_run_time(exp_dir_list, plot_against_wall_time=True, average_time_parameter_sweep=False, file_to_use_for_timing = 'logfile.000000.out'): +def calculate_month_run_time(exp_dir_list, loc_list, plot_against_wall_time=True, average_time_parameter_sweep=False, file_to_use_for_timing = 'logfile.000000.out'): """A script that takes a list of experiment names as input, and plots the time taken to run each month in that experiment vs the wall time. """ try: @@ -15,10 +15,17 @@ def calculate_month_run_time(exp_dir_list, plot_against_wall_time=True, average_ print('Environment variables GFDL_DATA must be set') sys.exit(0) + mean_time_arr = [] for exp_dir in exp_dir_list: - exp_dir_full = GFDL_DATA+'/'+exp_dir+'/' + exp_idx = exp_dir_list.index(exp_dir) + location = loc_list[exp_idx] + + if location=='GFDL_DATA': + exp_dir_full = GFDL_DATA+'/'+exp_dir+'/' + else: + exp_dir_full = location+'/'+exp_dir+'/' #Finds all the months for particular experiment months_to_check=os.listdir(exp_dir_full) @@ -68,19 +75,26 @@ def calculate_month_run_time(exp_dir_list, plot_against_wall_time=True, average_ else: plt.plot(month_num_arr[:-1], delta_t_arr, label=exp_dir) plt.xlabel('Month number') - + mean_time_arr.append(np.mean(delta_t_arr)) + plt.legend() plt.ylabel('Wall time elapsed per month (minutes)') if average_time_parameter_sweep: ax = plt.gca() - ax.set_xscale('log') + # ax.set_xscale('log') + + return mean_time_arr if __name__=="__main__": - exp_dir_list = ['socrates_test_mk44_chunk_size_longer_'+str(i) for i in [1, 2, 4,8,16,32,64, 128,256, 512]] + #exp_dir_list = ['socrates_ga3_chunk_size_'+str(i) for i in [4,8,16,32,64]] + + exp_dir_list = ['soc_dsa_files_smooth_topo_mk1', 'soc_ga3_files_smooth_topo_mk1', 'soc_ga7_files_smooth_topo_mk1','soc_ga3_files_smooth_topo_mk1_with_output_topo', 'soc_ga3_files_smooth_topo_old_fft_mk2_long', 'soc_ga3_files_smooth_topo_fftw_mk1_fresh_compile_long'] + + loc_list = ['GFDL_DATA'] * 4 + ['/scratch/sit204/mounts/isca_data/']*2 - calculate_month_run_time(exp_dir_list, plot_against_wall_time=False, average_time_parameter_sweep=True, file_to_use_for_timing='git_hash_used.txt') + mean_time_arr = calculate_month_run_time(exp_dir_list, loc_list, plot_against_wall_time=False, average_time_parameter_sweep=False, file_to_use_for_timing='git_hash_used.txt') plt.show() diff --git a/src/path_names b/src/path_names index 929ebc18c..3832506df 100644 --- a/src/path_names +++ b/src/path_names @@ -56,6 +56,7 @@ shared/diag_manager/diag_manager.F90 shared/diag_manager/diag_output.F90 shared/diag_manager/diag_table.F90 shared/diag_manager/diag_util.F90 +shared/fft/fftw.F90 shared/fft/fft99.F90 shared/fft/fft.F90 shared/field_manager/field_manager.F90 diff --git a/src/shared/fft/fft.F90 b/src/shared/fft/fft.F90 index 0878b6aff..10f3f9420 100644 --- a/src/shared/fft/fft.F90 +++ b/src/shared/fft/fft.F90 @@ -68,7 +68,9 @@ module fft_mod use fft99_mod, only: fft991, set99 #endif #endif - +#ifdef FFTW3 +use fftw3, only : fftw3_init, fftw3_end, grid_to_fourier_fftw, fourier_to_grid_fftw +#endif implicit none private @@ -188,6 +190,10 @@ module fft_mod !---------------------- private data ----------------------------------- +#ifdef FFTW3 + + +#endif ! tables for trigonometric constants and factors ! (not all will be used) real(R8_KIND), allocatable, dimension(:) :: table8 @@ -195,6 +201,9 @@ module fft_mod real , allocatable, dimension(:) :: table99 integer , allocatable, dimension(:) :: ifax + + + logical :: do_log =.true. integer :: leng, leng1, leng2, lenc ! related to transform size @@ -262,6 +271,7 @@ function fft_grid_to_fourier_float_2d (grid) result (fourier) ! argument "grid". ! !----------------------------------------------------------------------- + #ifdef SGICRAY # ifdef _CRAY ! local storage for cray fft @@ -303,10 +313,12 @@ function fft_grid_to_fourier_float_2d (grid) result (fourier) 'length of input data too small.') #endif !----------------------------------------------------------------------- -!----------------transform to fourier coefficients (+1)----------------- +!----------------transform to fourier coefficients (+1)----------------- num = size(grid,2) ! number of transforms - +#ifdef FFTW3 + call grid_to_fourier_fftw(num, leng1, lenc, grid, fourier) +#else #ifdef SGICRAY ! Cray/SGI fft scale = 1./real(leng) @@ -346,6 +358,7 @@ function fft_grid_to_fourier_float_2d (grid) result (fourier) enddo # endif #endif +#endif !----------------------------------------------------------------------- end function fft_grid_to_fourier_float_2d @@ -424,7 +437,9 @@ function fft_fourier_to_grid_float_2d (fourier) result (grid) #endif !----------------------------------------------------------------------- !----------------inverse transform to real space (-1)------------------- - +#ifdef FFTW3 + call fourier_to_grid_fftw(num, leng+1, lenc, fourier, grid) +#else #ifdef SGICRAY ! Cray/SGI fft scale = 1.0 @@ -469,6 +484,7 @@ function fft_fourier_to_grid_float_2d (fourier) result (grid) enddo # endif #endif +#endif !----------------------------------------------------------------------- @@ -502,6 +518,9 @@ function fft_grid_to_fourier_double_2d (grid) result (fourier) ! argument "grid". ! !----------------------------------------------------------------------- +#ifdef FFTW3 +! use fftw3 +#else #ifdef SGICRAY # ifdef _CRAY ! local storage for cray fft @@ -517,14 +536,15 @@ function fft_grid_to_fourier_double_2d (grid) result (fourier) # else ! local storage for temperton fft real, dimension(leng2,size(grid,2)) :: data - real, dimension(leng1,size(grid,2)) :: work -# endif + real, dimension(leng1,size(grid,2)) :: work +# endif #endif +#endif #if defined(SGICRAY) || defined(NAGFFT) real(R8_KIND) :: scale -#endif - integer :: j, k, num, len_grid +#endif + integer(kind=4) :: j, k, num, len_grid, i #ifdef NAGFFT integer :: ifail #endif @@ -549,9 +569,13 @@ function fft_grid_to_fourier_double_2d (grid) result (fourier) !----------------transform to fourier coefficients (+1)----------------- num = size(grid,2) ! number of transforms + +#ifdef FFTW3 +call grid_to_fourier_fftw(num, leng1, lenc, grid, fourier) +#else #ifdef SGICRAY ! Cray/SGI fft - scale = 1./float(leng) + scale = 1. / float(leng) # ifdef _CRAY call scfftm (-1,leng,num,scale, grid,leng1, fourier,lenc, & table8, work, 0) @@ -586,6 +610,7 @@ function fft_grid_to_fourier_double_2d (grid) result (fourier) enddo # endif #endif +#endif !----------------------------------------------------------------------- end function fft_grid_to_fourier_double_2d @@ -620,6 +645,9 @@ function fft_fourier_to_grid_double_2d (fourier) result (grid) ! argument "fourier". ! !----------------------------------------------------------------------- +#ifdef FFTW3 +! use fftw3 +#else #ifdef SGICRAY # ifdef _CRAY ! local storage for cray fft @@ -638,11 +666,12 @@ function fft_fourier_to_grid_double_2d (fourier) result (grid) real, dimension(leng1,size(fourier,2)) :: work # endif #endif +#endif #if defined(SGICRAY) || defined(NAGFFT) real(R8_KIND) :: scale #endif - integer :: j, k, num, len_fourier + integer(kind=4) :: j, k, num, len_fourier #ifdef NAGFFT integer :: ifail #endif @@ -667,7 +696,9 @@ function fft_fourier_to_grid_double_2d (fourier) result (grid) #endif !----------------------------------------------------------------------- !----------------inverse transform to real space (-1)------------------- - +#ifdef FFTW3 +call fourier_to_grid_fftw(num, leng+1, lenc, fourier, grid) +#else #ifdef SGICRAY ! Cray/SGI fft scale = 1.0 @@ -712,6 +743,7 @@ function fft_fourier_to_grid_double_2d (fourier) result (grid) enddo # endif #endif +#endif !----------------------------------------------------------------------- @@ -874,10 +906,14 @@ subroutine fft_init (n) ! variables that save length of transform leng = n; leng1 = n+1; leng2 = n+2; lenc = n/2+1 +#ifdef FFTW3 +call fftw3_init(leng, lenc) + +#else #ifdef SGICRAY # ifdef _CRAY ! initialization for cray -! float kind may not apply for cray +! float kind may not apply for cray allocate (table4(100+2*leng), table8(100+2*leng)) ! size may be too large? call scfftm (0,leng,1,0.0, dummy4, 1, cdummy4, 1, table4, dummy4, 0) call scfftm (0,leng,1,0.0, dummy8, 1, cdummy8, 1, table8, dummy8, 0) @@ -892,7 +928,7 @@ subroutine fft_init (n) # ifdef NAGFFT ! initialization for nag fft ifail8 = 0 - allocate (table8(100+2*leng)) ! size may be too large? + allocate (table8(100+2*leng)) ! size may be too large?a call c06fpf ( 1, leng, data8, 'i', table8, work8, ifail8 ) ! will not allow float kind for nag @@ -909,6 +945,7 @@ subroutine fft_init (n) allocate (ifax(10)) call set99 ( table99, ifax, leng ) # endif +#endif #endif module_is_initialized = .true. @@ -950,10 +987,12 @@ subroutine fft_end leng = 0; leng1 = 0; leng2 = 0; lenc = 0 + if (allocated(table4)) deallocate (table4) if (allocated(table8)) deallocate (table8) if (allocated(table99)) deallocate (table99) + module_is_initialized = .false. !----------------------------------------------------------------------- @@ -996,12 +1035,12 @@ program test allocate (ain(n+1,lot),aout(n+1,lot),four(n/2+1,lot)) call random_number (ain(1:n,:)) aout(1:n,:) = ain(1:n,:) - + call fft_init (n) ! transform grid to fourier and back four = fft_grid_to_fourier (aout) aout = fft_fourier_to_grid (four) - + ! print original and transformed do j=1,lot do i=1,n diff --git a/src/shared/fft/fftw.F90 b/src/shared/fft/fftw.F90 new file mode 100644 index 000000000..9a7619adf --- /dev/null +++ b/src/shared/fft/fftw.F90 @@ -0,0 +1,218 @@ +module fftw3 + use, intrinsic :: iso_c_binding + include 'fftw3.f03' + public :: fftw3_init + public :: fftw3_end + public :: grid_to_fourier_fftw + public :: fourier_to_grid_fftw + + interface grid_to_fourier_fftw + module procedure grid_to_fourier_double_2d_fftw + module procedure grid_to_fourier_float_2d_fftw + end interface + + interface fourier_to_grid_fftw + module procedure fourier_to_grid_double_2d_fftw + module procedure fourier_to_grid_float_2d_fftw + end interface + + private +# ifdef OVERLOAD_C8 +# define C_KIND_REAL C_DOUBLE +# define C_KIND_COMPLEX C_DOUBLE_COMPLEX +# endif + +# ifdef OVERLOAD_C4 +# define C_KIND_REAL C_FLOAT +# define C_KIND_COMPLEX C_FLOAT_COMPLEX +# endif + + + ! real -> complex + real(C_KIND_REAL), pointer :: real_input(:) + real(C_KIND_REAL), pointer :: real_output(:) + type(C_PTR) :: real_input_pointer, real_output_pointer + + ! complex -> real + complex(C_KIND_COMPLEX), pointer :: complex_input(:) + complex(C_KIND_COMPLEX), pointer :: complex_output(:) + type(C_PTR) :: complex_input_pointer, complex_output_pointer + + logical :: module_is_initialized = .false. +contains + +! +! Real to complex double precision fft +! +subroutine grid_to_fourier_double_2d_fftw(num, leng, lenc, grid, fourier) + +integer(kind=4), intent(in) :: num ! number of transformations +integer(kind=4), intent(in) :: leng ! length of each transformation +real(C_DOUBLE), intent(in) :: grid(leng, num) ! input grid +complex(C_DOUBLE_COMPLEX), intent(out) :: fourier(lenc, num) ! output fourier domain +real :: fact ! factor by which to scale results +integer :: i, j + +fact = 1.0 / (leng - 1) + +do j = 1, num + do i = 1, leng - 1 + real_input(i) = grid(i,j) + enddo +# ifdef OVERLOAD_C8 + call dfftw_execute_dft_r2c(real_input_pointer, real_input, complex_output) +# endif + do i = 1, lenc + fourier(i, j) = complex_output(i) * fact + enddo +enddo +return +end subroutine grid_to_fourier_double_2d_fftw + + +! +! Complex to real double precision fft +! +subroutine fourier_to_grid_double_2d_fftw(num, leng, lenc, fourier, grid) + +integer(kind=4), intent(in) :: num ! number of transformations +integer(kind=4), intent(in) :: leng ! length of each transformation +complex(C_DOUBLE_COMPLEX), intent(in) :: fourier(lenc, num) ! output fourier domain +real(C_DOUBLE), intent(out) :: grid(leng, num) ! input grid + +integer :: i, j + +do j = 1, num + do i = 1, lenc + complex_input(i) = fourier(i, j) + enddo +# ifdef OVERLOAD_C8 + call dfftw_execute_dft_c2r(complex_input_pointer, complex_input, real_output) +# endif + do i = 1, leng - 1 + grid(i, j) = real_output(i) + enddo + +enddo +return + +end subroutine fourier_to_grid_double_2d_fftw + +! +! Complex to real dingle precison fft +! +subroutine fourier_to_grid_float_2d_fftw(num, leng, lenc, fourier, grid) + +integer(kind=4), intent(in) :: num ! number of transformations +integer(kind=4), intent(in) :: leng ! length of each transformation +complex(C_FLOAT_COMPLEX), intent(in) :: fourier(lenc, num) ! output fourier domain +real(C_FLOAT), intent(out) :: grid(leng, num) ! input grid + +integer :: i, j + +do j = 1, num + do i = 1, lenc + complex_input(i) = fourier(i, j) + enddo +# ifdef OVERLOAD_C4 + call fftwf_execute_dft_c2r(complex_input_pointer, complex_input, real_output) +# endif + do i = 1, leng - 1 + grid(i, j) = real_output(i) + enddo + +enddo +return + +end subroutine fourier_to_grid_float_2d_fftw + +! +! Real to complex single precision fft +! +subroutine grid_to_fourier_float_2d_fftw(num, leng, lenc, grid, fourier) + +integer(kind=4), intent(in) :: num ! number of transformations +integer(kind=4), intent(in) :: leng ! length of each transformation +real(C_FLOAT), intent(in) :: grid(leng, num) ! input grid +complex(C_FLOAT_COMPLEX), intent(out) :: fourier(lenc, num) ! output fourier domain +real :: fact ! factor by which to scale results +integer :: i, j + +fact = 1.0 / (leng - 1) + +do j = 1, num + do i = 1, leng - 1 + real_input(i) = grid(i,j) + enddo +# ifdef OVERLOAD_C4 + call fftwf_execute_dft_r2c(real_input_pointer, real_input, complex_output) +# endif + do i = 1, lenc + fourier(i, j) = complex_output(i) * fact + enddo +enddo +return +end subroutine grid_to_fourier_float_2d_fftw + +! +! Initalise variables for fftw and setup plans +! +subroutine fftw3_init(leng, lenc) + integer, intent(in) :: leng, lenc + print *, 'Using FFTW3' + +# ifdef OVERLOAD_C8 + ! setup real -> complex transform + real_input_pointer = fftw_alloc_real(int(leng, C_SIZE_T)) + call c_f_pointer(real_input_pointer, real_input, [leng]) + + complex_output_pointer = fftw_alloc_complex(int(lenc, C_SIZE_T)) + call c_f_pointer(complex_output_pointer, complex_output, [leng]) + + call dfftw_plan_dft_r2c_1d(real_input_pointer, leng, real_input, complex_output, FFTW_ESTIMATE) + + ! setup complex -> real transform + complex_input_pointer = fftw_alloc_complex(int((lenc), C_SIZE_T)) + call c_f_pointer(complex_input_pointer, complex_input, [lenc]) + + real_output_pointer = fftw_alloc_real(int(leng, C_SIZE_T)) + call c_f_pointer(real_output_pointer, real_output, [leng]) + + call dfftw_plan_dft_c2r_1d(complex_input_pointer, leng, complex_input, real_output, FFTW_ESTIMATE) + +# else +# ifdef OVERLOAD_C4 + ! setup real -> complex transform + real_input_pointer = fftwf_alloc_real(int(leng, C_SIZE_T)) + call c_f_pointer(real_input_pointer, real_input, [leng]) + + complex_output_pointer = fftwf_alloc_complex(int(lenc, C_SIZE_T)) + call c_f_pointer(complex_output_pointer, complex_output, [leng]) + + call sfftw_plan_dft_r2c_1d(real_input_pointer, leng, real_input, complex_output, FFTW_ESTIMATE) + + ! setup complex -> real transform + complex_input_pointer = fftwf_alloc_complex(int((lenc), C_SIZE_T)) + call c_f_pointer(complex_input_pointer, complex_input, [lenc]) + + real_output_pointer = fftwf_alloc_real(int(leng, C_SIZE_T)) + call c_f_pointer(real_output_pointer, real_output, [leng]) + + call sfftw_plan_dft_c2r_1d(complex_input_pointer, leng, complex_input, real_output, FFTW_ESTIMATE) +# endif +# endif + + + module_is_initialized = .true. +end subroutine fftw3_init + +! +! Clean up all fftw variables +! +subroutine fftw3_end() + call fftw_cleanup() + module_is_initialized = .false. +end subroutine fftw3_end + + +end module fftw3 diff --git a/src/shared/mosaic/create_xgrid.c b/src/shared/mosaic/create_xgrid.c index e539a6ac1..50040f820 100644 --- a/src/shared/mosaic/create_xgrid.c +++ b/src/shared/mosaic/create_xgrid.c @@ -38,6 +38,7 @@ double dist_between_boxes(const double *x1, const double *y1, const double *z1, int inside_edge(double x0, double y0, double x1, double y1, double x, double y); int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double *intersect, double *u_a, double *u_q); +int isHeadNode(struct Node *list, struct Node nodeIn); /******************************************************************************* int get_maxxgrid @@ -1826,7 +1827,7 @@ double grid_box_radius(const double *x, const double *y, const double *z, int n) double dist_between_boxes(const double *x1, const double *y1, const double *z1, int n1, const double *x2, const double *y2, const double *z2, int n2) { - double dist; + double dist = 0; int i, j; for(i=0; i