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 startfi.nc file to apply fixed (= isothermal) surface temperatures across the planet.

Two options are available to do that:

Using newstart.e routine

       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

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:

rico_topography.png

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.

Using Python scripts

We have developed Python routines to easily modify start.nc and startfi.nc files. This requires the use of the xarray library.

• 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.

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('surface_mars.nc',decode_times=False) # can be any netcdf file (e.g. start/startfi.nc files)

# 2. WE READ THE VARIABLES
lat=nc['latitude']
lon=nc['longitude']
albedo=nc['albedo']
thermal_inertia=nc['thermal']
elevation=nc['zMOL']

# BELOW THE VARIABLES WE WANT TO UPDATE
new_elevation = np.empty((len(lat),len(lon)))
new_thermal_inertia = np.empty((len(lat),len(lon)))
new_albedo = np.empty((len(lat),len(lon)))

# LOOP TO MODIFY THE VARIABLES

# EXAMPLE 1 - CUSTOM TOPOGRAPHY
rico_topography=np.loadtxt('rico_topography.txt') # 360 lon x 180 lat
for i in range(0,len(lat),1):
    for j in range(0,len(lon),1):
        if(rico_topography[i,j]>1.):
            new_elevation[i,j]=rico_topography[i,j] # here you put whatever you want ; in this example, we use the topography data of rico_topography.txt
            new_albedo[i,j]=0.5 # here you put whatever you want ; we put high albedo because RICO is so bright his albedo must be high
            new_thermal_inertia[i,j]=2000. # here you put whatever you want
        else:
            new_elevation[i,j]=0. # here you put whatever you want
            new_albedo[i,j]=0.2 # here you put whatever you want
            new_thermal_inertia[i,j]=500. # here you put whatever you want
            
# EXAMPLE 2 - MARS WITH AN OCEAN
"""for i in range(0,len(lat),1):
    for j in range(0,len(lon),1):
        new_elevation[i,j]=max(elevation[i,j],-3.9) # here you put whatever you want ; in this example, we take the (MOLA) topography of present-day Mars and modify it to fill the topographic depressions with an ocean at an altitude of -3.9km
        new_albedo[i,j]=albedo[i,j] # here you put whatever you want
        new_thermal_inertia[i,j]=thermal_inertia[i,j] # here you put whatever you want"""

# SANITY CHECK PLOTS

fig1 = mpl.figure(1)
mpl.contourf(lon,lat,elevation)
mpl.xlabel('longitude (deg)')
mpl.ylabel('latitude (deg)')
mpl.show()

fig2 = mpl.figure(2)
mpl.contourf(lon,lat,new_elevation)
mpl.xlabel('longitude (deg)')
mpl.ylabel('latitude (deg)')
mpl.show()

# CREATE THE NEW NETCDF TOPOGRAPHY FILE

nc['albedo'].values = new_albedo
nc['zMOL'].values = new_elevation
nc['thermal'].values = new_thermal_inertia
nc.to_netcdf('new_topography.nc')