Modify start Files

From Planets
Revision as of 15:34, 7 November 2023 by Martin.turbet (talk | contribs)

Jump to: navigation, search

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 must modify the start.nc and startfi.nc files.

Two options are available to do that:

Using newstart.e routine

A few options to do that:

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 convert surface data files (elevation, albedo, thermal inertia) into GCM-compatible surface topography files. This requires the use of the xarray library.

Please find below examples: 

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')

Using a new topography in the model

Now that you have your netcdf surface topography file ready, simply add the surface topography netcdf file in your datagcm/surface_data/ directory, and then follow the tutorial here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Advanced_Use_of_the_GCM#How_to_Change_the_Topography_.28or_remove_it.29

Retrieved from "http://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php?title=Modify_start_Files&oldid=1744"