Modify start Files
Here we present a simple tutorial to change the initial conditions used in the Generic PCM. The initial conditions are stored in the start.nc and startfi.nc netCDF files, which are read by the Generic PCM. To modify the initial conditions, you thus must modify the start.nc and startfi.nc files.
In the examples below, we modify the start.nc and startfi.nc file to apply fixed (= isothermal) temperatures across the planet.
Two options are available to do that:
Using Python scripts
We have developed simple Python routines to easily modify start.nc and startfi.nc files. This requires the use of the xarray library.
In the example below, we modify the surface temperatures (stored in the startfi.nc file) to arbitrarily fix them to 350K in the North hemisphere, and 250K in the South hemisphere of the planet.
from numpy import *
import numpy as np
import matplotlib.pyplot as mpl
import math
import xarray as xr
# 1. WE GET THE SURFACE TOPOGRAPHY DATA
nc = xr.open_dataset('startfi.nc',decode_times=False)
# 2. WE READ THE VARIABLES
physical_points=nc['physical_points']
lat=nc['latitude']
lon=nc['longitude']
tsurf=nc['tsurf']
# BELOW THE VARIABLES WE WANT TO UPDATE
new_tsurf = np.empty(len(physical_points))
# LOOP TO MODIFY THE VARIABLES
# EXAMPLE - FIXED TEMPERATURES IN THE NORTH AND SOUTH HEMISPHERES
rico_topography=np.loadtxt('rico_topography.txt') # 360 lon x 180 lat
for i in range(0,len(physical_points),1):
if(lat[i]*180./np.pi >= 0.):
new_tsurf[I]=350. # north hemisphere
else:
new_tsurf[I]=250. # south hemisphere
# CREATE THE NEW NETCDF TOPOGRAPHY FILE
nc['tsurf'].values = new_tsurf
nc.to_netcdf('restartfi.nc')
The routine creates a restartfi.nc that can be copied into a startfi.nc and used as a new initial condition for the model.
However, to make more advanced modifications, it is sometimes necessary to switch to the other method (see below). For instance, to modify the 3D atmospheric temperature field, it is necessary to compute the potential temperature, which requires the use of the Exner function. As this operation is not trivial, we recommend to use the newstart.e routine (see below).
Using newstart.e routine
Newstart is an interactive tool to modify the start files (start.nc and startfi.nc). See https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Tool_Box#newstart:_a_fortran_program_to_modify_start_files for a tutorial on how to generate the executable file.
In the example below, we show how to fix all temperatures (atmosphere, surface, subsurface) to isothermal conditions, to a fixed temperature of 300K.
Once newstart.e is generated, you can execute the routine:
./newstart_XXX.e
and then type
0 # dd
When you execute newstart, you can use both a start2archive file or the start files (start.nc and startfi.nc). Then the interactive interface will propose to modify several physical quantities such as the gravity, the surface pressure or the rotation of the planet. At the end of the procedure, two files are created: restart.nc and restartfi.nc. They can be renamed and used as start files to initialize a new simulation.
- Create file start_archive.nc with start2archive.e compiled at grid resolution 64×48×32 using old file z2sig.def used previously
- Create files restart.nc and restartfi.nc with newstart.e compiled at grid resolution 32×24×25, using a new file z2sig.def (more details below on the choice of the z2sig.def).
- While executing newstart.e, you need to choose the answer '0 - from a file start_archive' and then press enter to all other requests.
else if (trim(modif) .eq. 'isotherm') then
write(*,*)'Isothermal temperature of the atmosphere,
& surface and subsurface'
write(*,*) 'Value of this temperature ? :'
203 read(*,*,iostat=ierr) Tiso
if(ierr.ne.0) goto 203
tsurf(1:ngridmx)=Tiso
tsoil(1:ngridmx,1:nsoilmx)=Tiso
Tset(1:iip1,1:jjp1,1:llm)=Tiso
flagtset=.true.
t(1:iip1,1:jjp1,1:llm)=Tiso
!! otherwise hydrost. integrations below
!! use the wrong temperature
!! -- NB: Tset might be useless
ucov(1:iip1,1:jjp1,1:llm)=0
vcov(1:iip1,1:jjm,1:llm)=0
q2(1:ngridmx,1:llm+1)=0
if (flagtset) then
DO l=1,llm
DO j=1,jjp1
DO i=1,iim
teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
ENDDO
teta (iip1,j,l)= teta (1,j,l)
ENDDO
ENDDO
A few options to do that:
- use existing topographic data for solar system objects. We have already formatted the data for most solar system objects, so please ask someone in the team if your topography map is not already here: https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/surface_data/
- use your own data files
You can also derive data files by digitalizing images. See below an example:
from numpy import *
import numpy as np
import matplotlib.pyplot as mpl
import math
from PIL import Image
# Here we first digitalize a png image name "rico_topography.png"
image = Image.open('rico_topography.png') # For simplicity, we work here in a 360 longitude pixels x 180 latitude pixels
pix = image.load()
image_width, image_height = image.size
topography_data=np.zeros((image_width, image_height),dtype='f')
# Here we convert the pixels into some rules for the elevation map
for i in range(0,image_width,1):
for j in range(0,image_height,1):
if(pix[i,j][0]<1. and pix[i,j][1]<1. and pix[i,j][2]<1.):
topography_data[i,j]=5.
else:
topography_data[i,j]=0.
# Here we write the topography map into an ascii file
name_file='rico_topography.txt'
with open(name_file, 'w') as f1:
for j in range(0,image_height,1):
for i in range(0,image_width-1,1):
print("{:12.5e}".format(topography_data[i,j]),file=f1,end=' ')
print("{:12.5e}".format(topography_data[image_width-1,j]),file=f1,end='\n')
This requires the use of the PIL library, and the rico_topography.png file.
