Isca
Isca copied to clipboard
ExeClim
Column model pr
This pull request adds single-column model functionality to the Isca 'framework' for modelling planetary atmospheres. When the column model is used, the dynamical core is bypassed but Isca's 'physics' routines can still be called.
This functionality is useful for testing new convection / radiation parametrizations as the column model runs a gazillion times faster (especially if you're only simulating one column) than the full GCM. Due to the way the code has been rewritten, the user can in principle simulate many columns (in lat and lon) at the same time, however currently only using one core. Attempting to use >1 core will cause the model to crash.
This development should very much be considered to be in the BETA phase, as it has not yet been sufficiently tested, and because it may in time become clear that some things implemented by this pull request should be done in a different way.
The column model is currently initiated as a bit of a hack. The line
from isca import ColumnCodeBase
in a python experiment configuration script sets a compiler flag -DCOLUMN_MODEL that tells the model to use the following files:
atmos_column/column.F90 atmos_column/column_grid.F90 atmos_column/column_init_cond.F90 atmos_column/column_initialize_fields.F90
to initialize the model (including constructing the model grid), do the model timestepping (using a leapfrog scheme as before), and handle input/output.
Some known issues at present:
- As mentioned previously, Isca running in column model 'mode' can only use one core.
- When configured as a single column,
interpolator_modwill not interpolate input netcdf files as intended, but instead return an array full of zeros. - Spatially varying topography cannot be specified, however the user can impose a spatially homogeneous non-zero geopotential height if they really want to.
Issue 1) likely arises because the code is trying to make use of the spec_mpp routines, in spite of the fact that when the column model is configured a spectral space is not defined. This issue in principle could be resolved by implementing the mpp rountines from an FMS gridpoint GCM, however I have yet to investigate this possibility sufficiently.
Issue 2) can probably be resolved by creating an interpolator_1D interface, however even in this scenario it is not clear to me what the desired functionality here would be (e.g. averaging over an input file with lat and lon? or just selecting the input data that corresponds to the single column's location in space?).
A hacky work around for Issue 2 is implemented in rrtm_radiation.f90 to enable the user to supply a single vertical profile for ozone via the rrtm_radiation namelist. This code works in an identical way to that in vert_coordinate.F90 where the user is able to supply pressure levels via the namelist.
Example test cases column_test.py and column_test_rrtm.py are provided to show how the model can be configured. Additionally supplied in support of column_test_rrtm.py is scm_interp_routine.py which can be used to interpolate model input files (e.g. ozone) to the model pressure levels, and select a single vertical profile (or create a globally averaged profile of the standard model input) before running the model (i.e. preprocessing). The example configured by column_test_rrtm.py makes use of this functionality to create a globally averaged vertical profile for ozone which is supplied to rrtm via the namelist.
It is worth mentioning the following:
The wind is prescribed (it needs to be non-zero at the surface to allow for latent and sensible surface heat fluxes). Currently the user can set a namelist variable 'surface_wind' that sets u_surf and v_surf = surface_wind / sqrt(2), so that wind_surf = sqrt(u_surf^2 + v_surf^2) = surface_wind. u and v at all other altitudes are set to zero (hardcoded).
At the moment the model needs to use the vertical turbulent diffusion parameterization in order for the mixed layer code to work. This is not very consistent as the u and v wind are prescribed and so the u,v tendenency from the diffusion is thrown away. Hence an implicit assumption when using the column model is that 'the dynamics' would restore the surface winds to their prescribed speed, so that du/dt total is zero. I have checked relatively extensively what using this parameterization is doing and it seems to be negligible / zero, so this is not really a concern. For example, when I tested the column by configuring it to solve for pure grey-gas radiative equilibirum, there was no appreciable difference between the result of the column model simulation and the analytical solution for grey gas radiative equilibrium.
Thanks for submitting this @ntlewis - there's a lot of great work here. I know we've already gone over this in person, I'll have another look at this and let you know what I think. :)
Is the code still essentially time-stepping to a solution? With each timestep then implementing a convective adjustment or similar (i.e. SBM or whatever)? If so, presumably the time to equilibration is set by the heat capacity the sound, and if that is small the model will equilibrate very quickly?
@gkvallis Yes the code is still time-stepping to a solution, using the same leapfrog method as used by the full GCM. On each timestep, convective adjustment, radiation and other physics routines are called as desired by the user (depending on the wishes of the user, the radiation may actually be called far less frequently, as in the full GCM). For an Earth-like atmosphere (i.e. 1 bar ish) you're correct that the model will equilibrate quickly if the surface heat capacity is small. I guess for a deeper atmosphere, say like Venus' it could take a little longer as radiative adjustment timescales in the deep atmosphere could be quite long. As the column model is computationally light, however, even long simulations take a very short time to run.
Hey @ntlewis - excited to use this for the first time.
I've just tried running the column_test.py for the first time, and it's not compiling. It gave me the following error. Do you think this is a problem that I'm having because this is my first time compiling it? The only thing I changed was I asked for 64 latitude columns instead of 1.
2019-01-17 16:57:02,775 - isca - INFO - mpifort -Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 -DCOLUMN_M
ODEL -I/usr/local/include -I/usr/local/include -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-intel -i4 -r8 -g -O2 -d
iag-disable 6843 -c /scratch/sit204/workdir_isca/codebase/_scratch_sit204_Isca_/code/src/atmos_column/column_initialize_fields.F90
2019-01-17 16:57:03,063 - isca - INFO - /scratch/sit204/workdir_isca/codebase/_scratch_sit204_Isca_/code/src/atmos_column/column_initi
alize_fields.F90(30): error #7002: Error in opening the compiled module file. Check INCLUDE paths. [TRANSFORMS_MOD]
2019-01-17 16:57:03,063 - isca - INFO - use transforms_mod, only: get_grid_domain
2019-01-17 16:57:03,063 - isca - INFO - ----------^
2019-01-17 16:57:03,064 - isca - INFO - /scratch/sit204/workdir_isca/codebase/_scratch_sit204_Isca_/code/src/atmos_column/column_initi
alize_fields.F90(30): error #6580: Name in only-list does not exist. [GET_GRID_DOMAIN]
2019-01-17 16:57:03,064 - isca - INFO - use transforms_mod, only: get_grid_domain
2019-01-17 16:57:03,064 - isca - INFO - --------------------------------^
2019-01-17 16:57:03,064 - isca - INFO - /scratch/sit204/workdir_isca/codebase/_scratch_sit204_Isca_/code/src/atmos_column/column_initi
alize_fields.F90(64): error #6406: Conflicting attributes or multiple declaration of name. [GET_GRID_DOMAIN]
2019-01-17 16:57:03,064 - isca - INFO - call get_grid_domain(is, ie, js, je)
2019-01-17 16:57:03,065 - isca - INFO - -----^
2019-01-17 16:57:03,068 - isca - INFO - compilation aborted for /scratch/sit204/workdir_isca/codebase/_scratch_sit204_Isca_/code/src/a
tmos_column/column_initialize_fields.F90 (code 1)
2019-01-17 16:57:03,074 - isca - INFO - make: *** [column_initialize_fields.o] Error 1
2019-01-17 16:57:03,074 - isca - INFO - ERROR: mkmf failed for column_isca.x
Exception in thread background thread for pid 11472:
Traceback (most recent call last):
File "/scratch/sit204/envs/pyi3/lib/python3.6/threading.py", line 916, in _bootstrap_inner
self.run()
File "/scratch/sit204/envs/pyi3/lib/python3.6/threading.py", line 864, in run
self._target(*self._args, **self._kwargs)
File "/scratch/sit204/envs/pyi3/lib/python3.6/site-packages/sh.py", line 1540, in wrap
fn(*args, **kwargs)
File "/scratch/sit204/envs/pyi3/lib/python3.6/site-packages/sh.py", line 2459, in background_thread
handle_exit_code(exit_code)
File "/scratch/sit204/envs/pyi3/lib/python3.6/site-packages/sh.py", line 2157, in fn
return self.command.handle_command_exit_code(exit_code)
File "/scratch/sit204/envs/pyi3/lib/python3.6/site-packages/sh.py", line 815, in handle_command_exit_code
raise exc
sh.ErrorReturnCode_1:
RAN: /bin/bash /scratch/sit204/workdir_isca/codebase/_scratch_sit204_Isca_/build/column_isca/compile.sh
STDOUT:
loadmodules for emps-gv machines
Compiling postprocessing tools
loadmodules for emps-gv machines
/scratch/sit204/workdir_isca/codebase/_scratch_sit204_Isca_/build/column_isca/path_names
......................................................................................................................................
............................................ Makefile is ready.
mpifort -I/usr/local/include -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-intel -i4 -r8 -g -O2 -diag-disable 6843 -
c /scratch/sit204/workdir_isca/codebase/_scratch_sit204_Isca_/code/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/par
kind.f90
mpifort -I/usr/local/include -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-int... (12811 more, please see e.stdout)
STDERR:
@sit23 Hi Stephen, thats my fault. In column_initialise_fields the import get_grid_domain from transforms_mod needs to be changed to the same subroutine import but from spec_mpp_mod. I think I didn't catch this because my column code didn't recompile when I removed transforms_mod (and other things associated with the dynamical core) from path_names.
@sit23 My latest commit should fix this. Thanks for testing. Unfortunately you might get a few more of these, as it seems to have revealed that not all of my code recompiled when I changed the list of files to compile in path_names. If this proves too much of a headache, we can just re-add the extra files to path_names.
@sit23 Ok, all should work now. All instances where transforms_mod was being used have now been replaced with imports from spec_mpp_mod and column_grid_mod. There is however a PROBLEM. When I made this fix in RRTM using #ifdef statements, the compile was failing almost as if the -DCOMPILER_FLAG was not being passed to the file rrtm_radiation.f90. This is because it wasn't. I have since checked, and all of the files appended with .f90 do not receive these user defined compiler flags. My most recent commit adds a duplicate file rrtm_radiation.F90 which the column model's path_names file directs it to, and with this it compiles. Fortunately this is the only .f90 file I needed to communicate my custom compiler flag to. At minimum this commit will require rrtm_radiation.f90 > rrtm_radiation.F90, which will need the path_names files for all other models to recognize this change. It would perhaps also be sensible to consider renaming all .f90 files .F90, however maybe not in this commit. For the time being I have left both the .f90 and .F90 rrtm_radiation files in the code, however one will need to be deleted before merging this P/R. Is there a way to rename files with git so that it doesn't see them as entirely new code?
I have two follow up comments: .F90 turns on a C-like pre-processor which enables #define statements. Second is to change rrtm_radiation.f90 to rrtm_radiation.F90 I think I should use git mv rrtm_radiation.f90 rrtm_radiation.F90 assuming we deem this is appropriate. And then of course update the path_names files for all of the other models that use RRTM, e.g. isca, accordingly.
Hey @ntlewis - I ran column_test.py for many columns and it worked well. Only problem I noticed was that latb(1) wasn't assigned (in fortran indexing terms). This meant panoply struggled to plot the data, as the first latitude had the wrong boundary values. I've made a little tiny fix here, which you could add manually to your branch and add to the P/R?
https://github.com/ExeClim/Isca/commit/b62624dd5e7ba725eb2d2cb4f1b764274d372c42
I'll have a think about the pre-processor things.
@ntlewis Hi Neil, I find that we can add -cpp flag to gfortran compiler and -fpp flag to Intel fortran compiler ifort to enable preprocessing for the files ending with .f90, so the compliers should deal with #ifdef things correctly after adding these flags. In addition, the -cpp or -fpp flags have been added to the Makefile templates (e.g., src/extra/python/isca/templates/mkmf.template.gfort and src/extra/python/isca/templates/mkmf.template.ia64), so I'm not sure why they don't work for your -D flags...
@lqxyz @sit23 Ok, so it seems you're right about the flag. But, as you've discovered, it is already included in the Isca compilation. It seems like -D flags aren't being passed to the compiler when .f90 files are compiled. The question to answer then is: Why aren't the -D flags being handed to the compiler when a .f90 file is being compiled?
Possibly because of this?
https://github.com/ExeClim/Isca/blob/9cfba874024d5215a7f53c235c541b3dabc6e6d3/bin/mkmf#L72
You could probably test either by directly modifying mkmf, or by adding "-cpp" to $OTHERFLAGS or equivalent.
@jamesp Would the sensible thing be to just hand the same flags to the compiler for both .F90 and .f90? As @lqxyz as pointed out, the -cpp/-fpp for gfortran/ifort should allow the use of these extra compiler flags.
It’s a good question. That the compiler doesn’t do it by default makes me think there could be specific cases where it’s not a good idea, but only knowing enough to be dangerous, I’m not sure what they are!
Worth at least trying the change and seeing if it still works, but it does seem the more “proper” approach would be to rename the files F90.
On 21 Jan 2019, at 12:18, Neil Lewis [email protected] wrote:
@jamesp Would the sensible thing be to just hand the same flags to the compiler for both .F90 and .f90? As @lqxyz as pointed out, the -cpp/-fpp for gfortran/ifort should allow the use of these extra compiler flags.
— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.
Yeah ok, I'll change the mkmf quickly and see if it still compiles, however I'd be inclined to move everything to .F90 in a different P/R, and then perhaps to move rrtm_radiation.f90 to .F90 in this one. Thoughts @sit23 ?
Will changing from f90 to F90 work properly on a Mac? The default file system on Macs is case insensitive so there is no difference in the .f90 and .F90 files themselves (and one cannot have both in the same directory).
Having said this, there may be a difference in what a compiler does if F90 or f90 are explicitly mentioned in the script.
On 21 Jan 2019, at 12:26 pm, Neil Lewis [email protected] wrote:
Yeah ok, I'll change the mkmf quickly and see if it still compiles, however I'd be inclined to move everything to .F90 in a different P/R, and then perhaps to move rrtm_radiation.f90 to .F90 in this one. Thoughts @sit23 https://github.com/sit23 ?
— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ExeClim/Isca/pull/82#issuecomment-456057190, or mute the thread https://github.com/notifications/unsubscribe-auth/AG0cn7RTNrrec0lYkKKZsn6O9wgbqQt9ks5vFbIGgaJpZM4aBN7P.
@gkvallis My guess is that the compiler knows the difference, even though the file system is not case sensitive. I guess this could be checked however with a simple test program (@lqxyz sent me one that would serve the purpose).
@ntlewis I think you are right that the compiler will know about case differences. Worth checking though.
I think a problem would only arise if one file existed and the compiler tried explicitly to do something with the other by name; for example, read .F90 and write .f90. Here are some command-line tests:
$ touch foo
$ ls FOO # succeeds, probably using stat() or similar
$ test -f FOO & echo yes # ditto
$ find . -name FOO # fails, probably readdir() (then case-sensitive)
$ cat foo >> FOO # grows forever unless initially empty! open() or similar
For now I've just resolved the compiler issue my moving rrtm_radiation.f90 to rrtm_radiation.F90. @sit23 you mentioned you had an issue where you could run with many latitudes and two longitudes, but not one. Would you be able to point me to a branch where this is reproduce-able so I can try and search for the issue? :smile:
@sit23 A while ago, I noticed that we needed to move rrtm from .f90 to .F90 for the column model's ifdef statements to work properly. At the time I made an error and just did mv rrtm_radiation.f90 rrtm_radiation.F90 and then git add rrtm_radiation.F90. It seems that this has removed the version history for rrtm, which I presume is a bad thing. Do you know a way to get it back? The actual file that's there is completely fine, its just lost all of the version history pre-dating my mv from .f90 to .F90. @dennissergeev this is the sort of thing you might be able to help with also...
@sit23 @lqxyz have you encountered any issues using the column model since January that you have local fixes for? If so could you send them to me so that I can put them in this pr? Other than the .f90 > .F90 issue (see above), I think this is pretty much ready to go. The only other thing I need to do is make a SocratesColumnCodeBase to use Column + Socrates, but this shouldn't take too long. I'll also make a test case for this.
As far as I understand, we can still track the changes for a renamed file (see e.g. https://stackoverflow.com/a/2314757/5365232).
As far as I understand, we can still track the changes for a renamed file (see e.g. https://stackoverflow.com/a/2314757/5365232).
@dennissergeev I don't think this is the case for what's happened here. I think I deleted the file and then made a new one instead, which was a mistake. For example, if you go to my fork and then look at the column_model_pr branch, the file rrtm_radiation.F90 doesn't have any of the old commit messages attached to it.
Do you know of a way to manually attach a 'different' file's history to a new one, in the case that they're actually the same file, but with different names?
Looking at the 'files changed' tab, I'm a little confused Github has picked up that the pr does rrtm_radiation.f90 -> rrtm_radiation.F90, and only indicates that 46 lines are added (as opposed to ~1000 odd if it was counting every line of the file). However, if you then click on 'view file' for rrtm_radiation.F90, the history only shows my two commits involving it.
So I read the thread above, but I am not convinced renaming the extension is the right solution here - after all, all other fortran files seem to have a lower-case extension - better to keep it this way?
Regarding the issue of deleting/moving the file, it shouldn't be a problem because git mv is just a shorthand for mv && git add. And replacing the file should be picked up automatically by git. Have you tried looking up the file's history via the command line (git log)?