Single point modeling of CLMU

Junjie Yu, 2023-12-06, Manchester, UK

ENV: RL_urban; docker container: clmu; docker image: clmu-app:1

1 Get the shell script for running CESM/CLM

simulation location : New York City

New York City coordinates can be replaced with other city coordinates

40°42′46″N 74°00′22″W from https://en.wikipedia.org/wiki/New_York_City

New York City = {‘lat’: 40.42, ‘lon’: -74.00}

This case was forced with data from the previous simulation under SSP3-7.0 scenario. ref

This case requires users to manually download certain forcing data and make changes to the case configuration. Therefore, the user can choose to practice with a CTSM history simulation, which will automatically download all the required data for the user.

from utils.cesm import *

print(now_time())

CASESCRIPT = "/p/project/clm5.0/cime/scripts"
CASESCRIPT_local = "/home/zhonghua/p/project/clm5.0/cime/scripts"
CASEROOT_local = "/home/zhonghua/p/scratch/CESMDATAROOT/CaseOutputs"
DOUT_S_ROOT = "/home/zhonghua/p/scratch/CESMDATAROOT/Archive"

# user defined script and json file
# parameters in the json file and scrpits shoule be consistent with PTS_script.sh and config_newyork.json
scriptpath = "./scr/PTS_SSP370.sh"
configpath = "config_SSP370.json"


cesm = cesm_run(CASESCRIPT, CASESCRIPT_local, CASEROOT_local, DOUT_S_ROOT, configpath)
#CASESCRIPT, the path of the CESM2/CTSM scripts folder
#    e.g., CASESCRIPT = "/p/project/myclm/cime/scripts"
#CASESCRIPT_local, the path of the CESM2/CTSM scripts folder in the local machine
#    e.g., CASESCRIPT_local = "/home/zhonghua/p/project/clm5.0/cime/scripts"
#CASEROOT_local, the path of the CESM2/CTSM case folder in the local machine
#    e.g., CASEROOT_local = "/home/zhonghua/p/scratch/CESMDATAROOT/CaseOutputs"
#DOUT_S_ROOT, the path of the CESM2/CTSM output folder
#    e.g., DOUT_S_ROOT = "/p/scratch/CESMDATAROOT/Archive"
#JSONPATH, the path of the json file or the json file
#    e.g., JSONPATH = "config.json"
#    e.g., JSONPATH = {"case_lat": 0.9,
#                        "case_lon": 1.25,
#                        "case_name": "test",
#                        "fsurdat": "/p/project/myclm/cime/scripts/f19_g16.B1850.f19_g16.c090918.clm2.r.0.5x0.5_20090202.nc"}


# password = "***"
password = input("Password of server: ")

# reset case
# if your server need password to run sudo command, you need to add "sudo -S" before the command and add the password
cesm.reset_case(password=password)

# change the permission of the script
command_0 = "sudo -S docker exec myclm chmod +x /p/pyscripts/PTS_SSP370.sh"
run_command(command_0, password)

# run the script
with open("/home/zhonghua/p/pyscripts/PTS_SSP370.sh", "w") as f:
    f.write(cesm.create_case(scriptpath))
command = "sudo -S docker exec myclm sh -c '/p/pyscripts/PTS_SSP370.sh' "
run_command(command, password)

d0 = cesm.nc_view()

Current time: 2023-12-09 16:26:05

[sudo] password for zhonghua: [sudo] password for zhonghua:

2 Modify the data

cesm.config['case_start_year'] = "2055"
cesm.config['case_start_year'] = "2055"
with open("/home/zhonghua/p/pyscripts/PTS_SSP370.sh", "w") as f:
    f.write(cesm.create_case(scriptpath))
run_command(command, password)
d2055= cesm.nc_view()

d2055

<xarray.Dataset>
Dimensions:             (levgrnd: 25, levlak: 10, levdcmp: 1, time: 3,
                         hist_interval: 2, lon: 1, lat: 1, column: 14,
                         gridcell: 1, landunit: 5, pft: 28)
Coordinates:
  * levgrnd             (levgrnd) float32 0.01 0.04 0.09 ... 19.48 28.87 42.0
  * levlak              (levlak) float32 0.05 0.6 2.1 4.6 ... 25.6 34.33 44.78
  * levdcmp             (levdcmp) float32 1.0
  * time                (time) object 2055-07-01 00:00:00 ... 2055-09-01 00:0...
  * lon                 (lon) float32 286.2
  * lat                 (lat) float32 40.05
Dimensions without coordinates: hist_interval, column, gridcell, landunit, pft
Data variables: (12/54)
    mcdate              (time) int32 ...
    mcsec               (time) int32 ...
    mdcur               (time) int32 ...
    mscur               (time) int32 ...
    nstep               (time) int32 ...
    time_bounds         (time, hist_interval) object ...
    ...                  ...
    pfts1d_itype_veg    (pft) int32 ...
    pfts1d_itype_col    (pft) int32 ...
    pfts1d_itype_lunit  (pft) int32 ...
    pfts1d_active       (pft) int32 ...
    TREFMNAV            (time, landunit) float32 ...
    TREFMXAV            (time, landunit) float32 ...
Attributes: (12/38)
    title:                                     CLM History file information
    comment:                                   NOTE: None of the variables ar...
    Conventions:                               CF-1.0
    history:                                   created on 12/09/23 16:28:13
    source:                                    Community Land Model CLM4.0
    hostname:                                  clm5-docker
    ...                                        ...
    ctype_urban_shadewall:                     73
    ctype_urban_impervious_road:               74
    ctype_urban_pervious_road:                 75
    cft_c3_crop:                               1
    cft_c3_irrigated:                          2
    time_period_freq:                          month_1

xarray.Dataset

• Dimensions:
  • levgrnd: 25
  • levlak: 10
  • levdcmp: 1
  • time: 3
  • hist_interval: 2
  • lon: 1
  • lat: 1
  • column: 14
  • gridcell: 1
  • landunit: 5
  • pft: 28
• Coordinates: (6)
  • levgrnd
    (levgrnd)
    float32
    0.01 0.04 0.09 ... 19.48 28.87 42.0

    long_name :

    coordinate soil levels

    units :

    m

    array([1.000000e-02, 4.000000e-02, 9.000000e-02, 1.600000e-01, 2.600000e-01,
           4.000000e-01, 5.800000e-01, 8.000000e-01, 1.060000e+00, 1.360000e+00,
           1.700000e+00, 2.080000e+00, 2.500000e+00, 2.990000e+00, 3.580000e+00,
           4.270000e+00, 5.060000e+00, 5.950000e+00, 6.940000e+00, 8.030000e+00,
           9.795000e+00, 1.332777e+01, 1.948313e+01, 2.887072e+01, 4.199844e+01],
          dtype=float32)

  • levlak
    (levlak)
    float32
    0.05 0.6 2.1 ... 25.6 34.33 44.78

    long_name :

    coordinate lake levels

    units :

    m

    array([ 0.05 ,  0.6  ,  2.1  ,  4.6  ,  8.1  , 12.6  , 18.6  , 25.6  , 34.325,
           44.775], dtype=float32)

  • levdcmp
    (levdcmp)
    float32
    1.0

    long_name :

    coordinate soil levels

    units :

    m

    array([1.], dtype=float32)

  • time
    (time)
    object
    2055-07-01 00:00:00 ... 2055-09-...

    long_name :

    time

    bounds :

    time_bounds

    array([cftime.DatetimeNoLeap(2055, 7, 1, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2055, 8, 1, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2055, 9, 1, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • lon
    (lon)
    float32
    286.2

    long_name :

    coordinate longitude

    units :

    degrees_east

    array([286.25], dtype=float32)

  • lat
    (lat)
    float32
    40.05

    long_name :

    coordinate latitude

    units :

    degrees_north

    array([40.052357], dtype=float32)

• Data variables: (54)
  • mcdate
    (time)
    int32
    ...

    long_name :

    current date (YYYYMMDD)

    [3 values with dtype=int32]

  • mcsec
    (time)
    int32
    ...

    long_name :

    current seconds of current date

    units :

    s

    [3 values with dtype=int32]

  • mdcur
    (time)
    int32
    ...

    long_name :

    current day (from base day)

    [3 values with dtype=int32]

  • mscur
    (time)
    int32
    ...

    long_name :

    current seconds of current day

    [3 values with dtype=int32]

  • nstep
    (time)
    int32
    ...

    long_name :

    time step

    [3 values with dtype=int32]

  • time_bounds
    (time, hist_interval)
    object
    ...

    long_name :

    history time interval endpoints

    [6 values with dtype=object]

  • date_written
    (time)
    |S16
    ...

    [3 values with dtype=|S16]

  • time_written
    (time)
    |S16
    ...

    [3 values with dtype=|S16]

  • area
    (lat, lon)
    float32
    ...

    long_name :

    grid cell areas

    units :

    km^2

    [1 values with dtype=float32]

  • landfrac
    (lat, lon)
    float32
    ...

    long_name :

    land fraction

    [1 values with dtype=float32]

  • landmask
    (lat, lon)
    float64
    ...

    long_name :

    land/ocean mask (0.=ocean and 1.=land)

    [1 values with dtype=float64]

  • pftmask
    (lat, lon)
    float64
    ...

    long_name :

    pft real/fake mask (0.=fake and 1.=real)

    [1 values with dtype=float64]

  • nbedrock
    (lat, lon)
    float64
    ...

    long_name :

    index of shallowest bedrock layer

    [1 values with dtype=float64]

  • ZSOI
    (levgrnd, column)
    float32
    ...

    long_name :

    soil depth

    units :

    m

    [350 values with dtype=float32]

  • DZSOI
    (levgrnd, column)
    float32
    ...

    long_name :

    soil thickness

    units :

    m

    [350 values with dtype=float32]

  • WATSAT
    (levgrnd, column)
    float32
    ...

    long_name :

    saturated soil water content (porosity)

    units :

    mm3/mm3

    [350 values with dtype=float32]

  • SUCSAT
    (levgrnd, column)
    float32
    ...

    long_name :

    saturated soil matric potential

    units :

    mm

    [350 values with dtype=float32]

  • BSW
    (levgrnd, column)
    float32
    ...

    long_name :

    slope of soil water retention curve

    units :

    unitless

    [350 values with dtype=float32]

  • HKSAT
    (levgrnd, column)
    float32
    ...

    long_name :

    saturated hydraulic conductivity

    units :

    mm s-1

    [350 values with dtype=float32]

  • ZLAKE
    (levlak, column)
    float32
    ...

    long_name :

    lake layer node depth

    units :

    m

    [140 values with dtype=float32]

  • DZLAKE
    (levlak, column)
    float32
    ...

    long_name :

    lake layer thickness

    units :

    m

    [140 values with dtype=float32]

  • grid1d_lon
    (gridcell)
    float64
    ...

    long_name :

    gridcell longitude

    units :

    degrees_east

    [1 values with dtype=float64]

  • grid1d_lat
    (gridcell)
    float64
    ...

    long_name :

    gridcell latitude

    units :

    degrees_north

    [1 values with dtype=float64]

  • grid1d_ixy
    (gridcell)
    int32
    ...

    long_name :

    2d longitude index of corresponding gridcell

    [1 values with dtype=int32]

  • grid1d_jxy
    (gridcell)
    int32
    ...

    long_name :

    2d latitude index of corresponding gridcell

    [1 values with dtype=int32]

  • land1d_lon
    (landunit)
    float64
    ...

    long_name :

    landunit longitude

    units :

    degrees_east

    [5 values with dtype=float64]

  • land1d_lat
    (landunit)
    float64
    ...

    long_name :

    landunit latitude

    units :

    degrees_north

    [5 values with dtype=float64]

  • land1d_ixy
    (landunit)
    int32
    ...

    long_name :

    2d longitude index of corresponding landunit

    [5 values with dtype=int32]

  • land1d_jxy
    (landunit)
    int32
    ...

    long_name :

    2d latitude index of corresponding landunit

    [5 values with dtype=int32]

  • land1d_wtgcell
    (landunit)
    float64
    ...

    long_name :

    landunit weight relative to corresponding gridcell

    [5 values with dtype=float64]

  • land1d_ityplunit
    (landunit)
    int32
    ...

    long_name :

    landunit type (vegetated,urban,lake,wetland,glacier or glacier_mec)

    [5 values with dtype=int32]

  • land1d_active
    (landunit)
    int32
    ...

    long_name :

    true => do computations on this landunit

    flag_values :

    [0 1]

    flag_meanings :

    FALSE TRUE

    valid_range :

    [0 1]

    [5 values with dtype=int32]

  • cols1d_lon
    (column)
    float64
    ...

    long_name :

    column longitude

    units :

    degrees_east

    [14 values with dtype=float64]

  • cols1d_lat
    (column)
    float64
    ...

    long_name :

    column latitude

    units :

    degrees_north

    [14 values with dtype=float64]

  • cols1d_ixy
    (column)
    int32
    ...

    long_name :

    2d longitude index of corresponding column

    [14 values with dtype=int32]

  • cols1d_jxy
    (column)
    int32
    ...

    long_name :

    2d latitude index of corresponding column

    [14 values with dtype=int32]

  • cols1d_wtgcell
    (column)
    float64
    ...

    long_name :

    column weight relative to corresponding gridcell

    [14 values with dtype=float64]

  • cols1d_wtlunit
    (column)
    float64
    ...

    long_name :

    column weight relative to corresponding landunit

    [14 values with dtype=float64]

  • cols1d_itype_col
    (column)
    int32
    ...

    long_name :

    column type (see global attributes)

    [14 values with dtype=int32]

  • cols1d_itype_lunit
    (column)
    int32
    ...

    long_name :

    column landunit type (vegetated,urban,lake,wetland,glacier or glacier_mec)

    [14 values with dtype=int32]

  • cols1d_active
    (column)
    int32
    ...

    long_name :

    true => do computations on this column

    flag_values :

    [0 1]

    flag_meanings :

    FALSE TRUE

    valid_range :

    [0 1]

    [14 values with dtype=int32]

  • pfts1d_lon
    (pft)
    float64
    ...

    long_name :

    pft longitude

    units :

    degrees_east

    [28 values with dtype=float64]

  • pfts1d_lat
    (pft)
    float64
    ...

    long_name :

    pft latitude

    units :

    degrees_north

    [28 values with dtype=float64]

  • pfts1d_ixy
    (pft)
    int32
    ...

    long_name :

    2d longitude index of corresponding pft

    [28 values with dtype=int32]

  • pfts1d_jxy
    (pft)
    int32
    ...

    long_name :

    2d latitude index of corresponding pft

    [28 values with dtype=int32]

  • pfts1d_wtgcell
    (pft)
    float64
    ...

    long_name :

    pft weight relative to corresponding gridcell

    [28 values with dtype=float64]

  • pfts1d_wtlunit
    (pft)
    float64
    ...

    long_name :

    pft weight relative to corresponding landunit

    [28 values with dtype=float64]

  • pfts1d_wtcol
    (pft)
    float64
    ...

    long_name :

    pft weight relative to corresponding column

    [28 values with dtype=float64]

  • pfts1d_itype_veg
    (pft)
    int32
    ...

    long_name :

    pft vegetation type

    [28 values with dtype=int32]

  • pfts1d_itype_col
    (pft)
    int32
    ...

    long_name :

    pft column type (see global attributes)

    [28 values with dtype=int32]

  • pfts1d_itype_lunit
    (pft)
    int32
    ...

    long_name :

    pft landunit type (vegetated,urban,lake,wetland,glacier or glacier_mec)

    [28 values with dtype=int32]

  • pfts1d_active
    (pft)
    int32
    ...

    long_name :

    true => do computations on this pft

    flag_values :

    [0 1]

    flag_meanings :

    FALSE TRUE

    valid_range :

    [0 1]

    [28 values with dtype=int32]

  • TREFMNAV
    (time, landunit)
    float32
    ...

    long_name :

    daily minimum of average 2-m temperature

    units :

    K

    cell_methods :

    time: mean

    [15 values with dtype=float32]

  • TREFMXAV
    (time, landunit)
    float32
    ...

    long_name :

    daily maximum of average 2-m temperature

    units :

    K

    cell_methods :

    time: mean

    [15 values with dtype=float32]

• Indexes: (6)
  • levgrnd
    PandasIndex

    PandasIndex(Index([0.009999999776482582,  0.03999999910593033,  0.09000000357627869,
             0.1599999964237213,  0.25999999046325684,   0.4000000059604645,
             0.5799999833106995,    0.800000011920929,    1.059999942779541,
             1.3600000143051147,   1.7000000476837158,   2.0799999237060547,
                            2.5,    2.990000009536743,   3.5799999237060547,
              4.269999980926514,    5.059999942779541,    5.949999809265137,
              6.940000057220459,    8.029999732971191,    9.795000076293945,
             13.327767372131348,   19.483129501342773,   28.870723724365234,
             41.998435974121094],
          dtype='float32', name='levgrnd'))

  • levlak
    PandasIndex

    PandasIndex(Index([0.05000000074505806,  0.6000000238418579,  2.0999999046325684,
             4.599999904632568,   8.100000381469727,  12.600000381469727,
            18.600000381469727,  25.600000381469727,   34.32500076293945,
            44.775001525878906],
          dtype='float32', name='levlak'))

  • levdcmp
    PandasIndex

    PandasIndex(Index([1.0], dtype='float32', name='levdcmp'))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([2055-07-01 00:00:00, 2055-08-01 00:00:00, 2055-09-01 00:00:00],
                dtype='object', length=3, calendar='noleap', freq='MS'))

  • lon
    PandasIndex

    PandasIndex(Index([286.25], dtype='float32', name='lon'))

  • lat
    PandasIndex

    PandasIndex(Index([40.052357], dtype='float32', name='lat'))

• Attributes: (38)

  title :

  CLM History file information

  comment :

  NOTE: None of the variables are weighted by land fraction!

  Conventions :

  CF-1.0

  history :

  created on 12/09/23 16:28:13

  source :

  Community Land Model CLM4.0

  hostname :

  clm5-docker

  username :

  root

  version :

  release-clm5.0.35-6-gd010e00d5

  revision_id :

  $Id: histFileMod.F90 42903 2012-12-21 15:32:10Z muszala $

  case_title :

  UNSET

  case_id :

  SSP370_exp

  Surface_dataset :

  surfdata_0.9x1.25_hist_16pfts_Irrig_CMIP6_simyr1850_c190214.nc

  Initial_conditions_dataset :

  arbitrary initialization

  PFT_physiological_constants_dataset :

  clm5_params.c171117.nc

  ltype_vegetated_or_bare_soil :

  1

  ltype_crop :

  2

  ltype_UNUSED :

  3

  ltype_landice_multiple_elevation_classes :

  4

  ltype_deep_lake :

  5

  ltype_wetland :

  6

  ltype_urban_tbd :

  7

  ltype_urban_hd :

  8

  ltype_urban_md :

  9

  ctype_vegetated_or_bare_soil :

  1

  ctype_crop :

  2

  ctype_crop_noncompete :

  2*100+m, m=cft_lb,cft_ub

  ctype_landice :

  3

  ctype_landice_multiple_elevation_classes :

  4*100+m, m=1,glcnec

  ctype_deep_lake :

  5

  ctype_wetland :

  6

  ctype_urban_roof :

  71

  ctype_urban_sunwall :

  72

  ctype_urban_shadewall :

  73

  ctype_urban_impervious_road :

  74

  ctype_urban_pervious_road :

  75

  cft_c3_crop :

  1

  cft_c3_irrigated :

  2

  time_period_freq :

  month_1

3 View the results

import matplotlib.pyplot as plt

print(now_time())
fig, ax = plt.subplots(1, 1, figsize=(10, 5))


d0['TREFMXAV'].assign_coords({'landunit': d0['land1d_ityplunit'].values}).mean(dim='time').plot.scatter(s=500, ax=ax, label="TREFMXAV 2035")
d0['TREFMNAV'].assign_coords({'landunit': d0['land1d_ityplunit'].values}).mean(dim='time').plot.scatter(s=500, ax=ax, label="TREFMNAV 2035")

d2055['TREFMXAV'].assign_coords({'landunit': d0['land1d_ityplunit'].values}).mean(dim='time').plot.scatter(s=500, ax=ax, label="TREFMXAV 2055")
d2055['TREFMNAV'].assign_coords({'landunit': d0['land1d_ityplunit'].values}).mean(dim='time').plot.scatter(s=500, ax=ax, label="TREFMNAV 2055")
ax.set_xlabel('Lanuint ID')
ax.set_ylabel('2m air temperature [C]')
ax.legend(loc='upper left', bbox_to_anchor=(1.05, 1))

# Add grid lines
ax.grid(color='gray', linestyle='dashed')

# Set grid lines to align with x-axis ticks
ax.set_axisbelow(True)

plt.show()

Current time: 2023-12-09 16:29:16