Difference between revisions of "Modify start Files"

From Planets
Jump to: navigation, search
(→‎Using newstart.e routine)
Line 1: Line 1:
 
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.
 
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.
+
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:
 
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.
 +
 +
<syntaxhighlight lang="python">
 +
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')
 +
 +
</syntaxhighlight>
 +
 +
The routine creates a ''restartfi.nc'' that can be copied into a ''startfi.nc'' and used as a new initial condition for the model.
  
 
== Using newstart.e routine ==
 
== Using newstart.e routine ==
 +
 +
 +
* 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
 
         else if (trim(modif) .eq. 'isotherm') then
Line 86: Line 136:
  
 
This requires the use of the PIL library, and the rico_topography.png file.
 
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.
 
−
 
−
<syntaxhighlight lang="python">
 
−
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')
 
−
 
−
</syntaxhighlight>
 

Revision as of 16:21, 7 November 2023

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.

Using newstart.e routine

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

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.

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