Difference between revisions of "Modify start Files"

From Planets
Jump to: navigation, search
(Using newstart.e routine)
(Using newstart.e routine)
 
(15 intermediate revisions by the same user not shown)
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 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
 
 
 
      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:
 
 
[[File:rico_topography.png|thumb|rico_topography.png]]
 
 
<syntaxhighlight lang="python">
 
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')
 
</syntaxhighlight>
 
 
 
This requires the use of the PIL library, and the rico_topography.png file.
 
  
 
== Using Python scripts ==
 
== 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.
+
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 simple 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.
 
 
* 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">
 
<syntaxhighlight lang="python">
Line 105: Line 19:
  
 
# 1. WE GET THE SURFACE TOPOGRAPHY DATA
 
# 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)
+
nc = xr.open_dataset('startfi.nc',decode_times=False)
  
 
# 2. WE READ THE VARIABLES
 
# 2. WE READ THE VARIABLES
 +
physical_points=nc['physical_points']
 
lat=nc['latitude']
 
lat=nc['latitude']
 
lon=nc['longitude']
 
lon=nc['longitude']
albedo=nc['albedo']
+
tsurf=nc['tsurf']
thermal_inertia=nc['thermal']
 
elevation=nc['zMOL']
 
  
 
# BELOW THE VARIABLES WE WANT TO UPDATE
 
# BELOW THE VARIABLES WE WANT TO UPDATE
new_elevation = np.empty((len(lat),len(lon)))
+
new_tsurf = np.empty(len(physical_points))
new_thermal_inertia = np.empty((len(lat),len(lon)))
 
new_albedo = np.empty((len(lat),len(lon)))
 
  
 
# LOOP TO MODIFY THE VARIABLES
 
# LOOP TO MODIFY THE VARIABLES
 +
# EXAMPLE - FIXED TEMPERATURES IN THE NORTH AND SOUTH HEMISPHERES
 +
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
  
# EXAMPLE 1 - CUSTOM TOPOGRAPHY
+
# CREATE THE NEW NETCDF TOPOGRAPHY FILE
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
+
nc['tsurf'].values = new_tsurf
 +
nc.to_netcdf('restartfi.nc')
  
fig1 = mpl.figure(1)
+
</syntaxhighlight>
mpl.contourf(lon,lat,elevation)
 
mpl.xlabel('longitude (deg)')
 
mpl.ylabel('latitude (deg)')
 
mpl.show()
 
  
fig2 = mpl.figure(2)
+
The routine creates a ''restartfi.nc'' that can be copied into a ''startfi.nc'' and used as a new initial condition for the model.
mpl.contourf(lon,lat,new_elevation)
 
mpl.xlabel('longitude (deg)')
 
mpl.ylabel('latitude (deg)')
 
mpl.show()
 
  
# CREATE THE NEW NETCDF TOPOGRAPHY FILE
+
However, to make more advanced modifications, it is sometimes necessary to switch to another, more advanced method (see below). For instance, to modify the 3D atmospheric temperature field, it is necessary to compute and modify 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.
  
nc['albedo'].values = new_albedo
+
Once ''newstart.e'' is generated, you can execute the routine:
nc['zMOL'].values = new_elevation
+
<syntaxhighlight lang="bash">
nc['thermal'].values = new_thermal_inertia
+
./newstart_XXX.e # with XXX depending on the options used at the compilation stage
nc.to_netcdf('new_topography.nc')
+
</syntaxhighlight>
  
 +
and then type
 +
<syntaxhighlight lang="bash">
 +
1 # to modify start.nc and startfi.nc files
 +
n
 +
# press enter to pass
 +
isotherm # this selects the option to fix all temperatures to isothermal conditions
 +
300. # value of temperature
 +
# press enter to validate
 
</syntaxhighlight>
 
</syntaxhighlight>
 +
 +
This routine creates ''restart.nc'' and ''restartfi.nc'' files that can be copied into ''start.nc'' and ''startfi.nc'', respectively, and used as a new initial conditions for the Generic PCM.
 +
 +
Note 1: we invite you to go through all existing options of ''newstart.e'' to get an idea of which variables can be modified.
 +
Note 2: if the option you want does not exist, you can create your own in the ''newstart.e'' routine by modifying the LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F fortran file. You can easily copy and adapt existing examples. Once the changes have been made, you need to recompile the ''newstart.e'' routine and run it!

Latest revision as of 16:49, 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 simple 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
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 another, more advanced method (see below). For instance, to modify the 3D atmospheric temperature field, it is necessary to compute and modify 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 # with XXX depending on the options used at the compilation stage

and then type

1 # to modify start.nc and startfi.nc files
n
 # press enter to pass
isotherm # this selects the option to fix all temperatures to isothermal conditions
300. # value of temperature
 # press enter to validate

This routine creates restart.nc and restartfi.nc files that can be copied into start.nc and startfi.nc, respectively, and used as a new initial conditions for the Generic PCM.

Note 1: we invite you to go through all existing options of newstart.e to get an idea of which variables can be modified. Note 2: if the option you want does not exist, you can create your own in the newstart.e routine by modifying the LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F fortran file. You can easily copy and adapt existing examples. Once the changes have been made, you need to recompile the newstart.e routine and run it!