ocean.org

ModelE Land Ice (LIME) Ocean Coupling

Table of contents

• Introduction
• Summary
• ModelE ocean runs
  • Build rundeck and compile
    • E1oM20
  • First run
    • Run 1 hour
    • Run longer until {YEAR,MONTH,DATE,HOUR}E defined in I
  • Postprocess
• ModelE to PISM
  • Post-process one year data from ModelE to more generic format
  • Set up ModelE domain
  • Import
    • Meta
    • Temperature & salinity
  • Reproject
    • Greenland
    • Antarctica
  • Export to NetCDF
    • Display: Antarctic salt at all depths at time 0
    • Display: Greenland temp at all times at depth 0
    • Merge ModelE output with LIME metadata
    • Select appropriate T & S depth for each cell
    • NOTDONE Select appropriate T & S for each basin
    • Create basin IDs
    • Reformat to meet PISM requirements
  • Test run in PISM

Introduction

Notes, code, and documentation on coupling the ModelE ocean with PISM.

Eventual goal is bi-directional coupling so that ModelE ocean (temperature, salinity) forces PISM, and PISM freshwater is returned to ModelE.

At the moment coupling is one-way: We drive PISM with ModelE, using monthly coupling.

Summary

This summary may be out of date. See TOC, section headings, and code within for definitive code used to generate the products. Work plan is, roughly:

• Select ModelE ocean cells bordering Antarctica and Greenland
  • Extract pot_temp, salinity, and other needed metadata
• Reproject from ModelE grid to PISM grid
• Select relevant depth of ocean properties
  • GZ draft, bathymetric routing to shelf break, etc.

WARN

Use nearest neighbor sampling otherwise interpolation would occur, which might be OK for T or S, but unlikely to produce reasonable density which is f(T,S).

ModelE ocean runs

Template	Comment
E4M20	4x5 resolution qflux ocean
E1oM20	4x5 dynamic ocean
E6F40	2x2.5 resolution qflux ocean
E1F40o13	2x2.5 dynamic 13 layer ocean
E1F40o32	2x2.5 dynamic 32 layer ocean

Build rundeck and compile

E1oM20

cd ~/projects/GISS/modelE_2.1_branch/decks

RUNNAME=E1oM20_r01

make rundeck RUN=${RUNNAME} RUNSRC=E1oM20 OVERWRITE=YES

 # YEARI=1999,MONTHI=12,DATEI=1,HOURI=0,
 # YEARE=2011,MONTHE=1,DATEE=1,HOURE=0,     KDIAG=13*0,
 # ISTART=2,IRANDI=0, YEARE=1999,MONTHE=12,DATEE=1,HOURE=1,JWRITE=1
                                 
make clean RUN=${RUNNAME}
make -j setup RUN=${RUNNAME} EXTRA_FFLAGS="-Ddbgcsv_silent"

First run

Run 1 hour

../exec/runE ${RUNNAME} -cold-restart -np 2

Run longer until {YEAR,MONTH,DATE,HOUR}E defined in I

../exec/runE ${RUNNAME} -np 2

Examine outputs in run folder…

Postprocess

• From decks folder

cd ${RUNNAME}/

mkdir post
cd post

# annual sum of first 20 years
sumfiles ../???10{0,1}[0-9]* # ANN1000-1019.acc${RUNNAME}.nc

# monthly sum of first 20 years
MONTHS=$(locale mon|tr '[:lower:]' '[:upper:]'| tr ';' '\n' |cut -c1-3)
YEARS=$(seq 1000 1019|tr '\n' ','| sed 's/,$//')
parallel "sumfiles ../{1}10[0-1][0-9]*" ::: ${MONTHS}


scaleacc ANN1000-1019.acc${RUNNAME}.nc aij,aj,ajl,aijl,aijk,ijhc
ncview ANN1000-1019.aij${RUNNAME}.nc
# look at "impm_lndice", "MICB", etc.

ModelE to PISM

Post-process one year data from ModelE to more generic format

RUNNAME=E1oM20_r01
YEAR=2001

month_names=($(locale mon|tr '[:lower:]' '[:upper:]'| tr ';' '\n' |cut -c1-3))
for m in ${month_names[@]}; do
  scaleacc ${CMRUNDIR}/${RUNNAME}/${m}${YEAR}.acc${RUNNAME}.nc oijl
done

# just need one of these for metadata
scaleacc ${CMRUNDIR}/${RUNNAME}/JAN${YEAR}.acc${RUNNAME}.nc aij
mv JAN${YEAR}.aij${RUNNAME}.nc aij.nc

Set up ModelE domain

• Note, 4x5 arrays have a smaller polar cell.
• Easiest to set two domains
  • 4x5 +- 92 °
  • 2x2.5 +- 90 °

Resample to the +- 90 ° domain.

grass ./G_ModelE/PERMANENT
# g.region -ps n=90 s=-90 w=-180 e=180 nsres=1 ewres=1

g.mapset -c 4x5

# hack for ModelE polar cells
g.region -p n=92 s=-92 w=-180 e=180 nsres=4 ewres=5
g.region save=4x5_92

g.mapset -c 2x2p5
g.region -p n=90 s=-90 nsres=2 ewres=2.5
g.region save=2x2p5_90

Import

Meta

g.mapset 4x5
r.in.gdal -o input=NetCDF:aij.nc:landicefr output=landicefr
r.in.gdal -o input=NetCDF:aij.nc:landfr output=landfr
r.in.gdal -o input=NetCDF:aij.nc:ocnfr output=ocnfr

g.mapset 2x2p5
r.resample input=landicefr@4x5 output=landicefr
r.resample input=landfr@4x5 output=landfr
r.resample input=ocnfr@4x5 output=ocnfr

Temperature & salinity

RUNNAME=E1oM20_r01
YEAR=2001

nlevel=$(ncks -M -C -v zoc JAN${YEAR}.oijl${RUNNAME}.nc |grep "zoc = "|tr -cd '[:digit:]')
levels=$(seq -w ${nlevel})

month_names=$(locale mon|tr '[:lower:]' '[:upper:]'| tr ';' '\n' |cut -c1-3)
month_nums=$(seq -w 1 12)

g.mapset 4x5
parallel "r.in.gdal -o input=NetCDF:{1}${YEAR}.oijl${RUNNAME}.nc:{4} output={=4 s/.*_// =}_m{2}_z{3} band={3}" ::: ${month_names} :::+ ${month_nums} ::: ${levels} ::: salt pot_temp

g.mapset 2x2p5
parallel "r.resample input={3}_m{1}_z{2}@4x5 output={3}_m{1}_z{2}" ::: ${month_nums} ::: ${levels} ::: salt temp

Reproject

Greenland

The r.grow.distance flood-fills the values with nearest neighbor. This is because ModelE coast may not match PISM coast, but we will want to be able to access the T and S values nearest the coast.

grass ./G_GL/PERMANENT

RUNNAME=E1oM20_r01
YEAR=2001

nlevel=$(ncks -M -C -v zoc JAN${YEAR}.oijl${RUNNAME}.nc |grep "zoc = "|tr -cd '[:digit:]')
levels=$(seq -w ${nlevel})

month_names=$(locale mon|tr '[:lower:]' '[:upper:]'| tr ';' '\n' |cut -c1-3)
month_nums=$(seq -w 1 12)

r.proj location=G_ModelE mapset=2x2p5 input=landicefr method=nearest 
r.proj location=G_ModelE mapset=2x2p5 input=landfr method=nearest
r.proj location=G_ModelE mapset=2x2p5 input=ocnfr method=nearest

rasters=$(r.proj location=G_ModelE mapset=2x2p5 -l |grep -E '^temp_|^salt_')
parallel --progress --bar "r.proj -n --q location=G_ModelE mapset=2x2p5 input={1} output={1}_E method=nearest" ::: ${rasters}
parallel --progress --bar "r.grow.distance --q input={1}_E value={1}" ::: ${rasters}

# for var in temp salt; do
#   for z in ${levels}; do
#     rasters=$(g.list type=raster pattern="${var}_m??_*z${z}" sep=,)
#     r.to.rast3 --q input=${rasters} output=${var}_z${z}
#   done
# done

# Generate 12 3D rasters: For each month, 1 3D raster of all levels
g.region t=${nlevel} -p3
parallel "r.to.rast3 --q input=\$(g.list type=raster pattern=\"{1}_m{2}_z[0-9]?\" sep=,) output={1}_m{2}" ::: temp salt ::: ${month_nums}

# parallel "r3.out.netcdf -p input={1}_m{2} output={1}_m{2}_GL.nc" ::: temp salt ::: ${month_nums}

exit

Antarctica

grass ./G_AQ/PERMANENT

RUNNAME=E1oM20_r01
YEAR=2001

nlevel=$(ncks -M -C -v zoc JAN${YEAR}.oijl${RUNNAME}.nc |grep "zoc = "|tr -cd '[:digit:]')
levels=$(seq -w ${nlevel})

month_names=$(locale mon|tr '[:lower:]' '[:upper:]'| tr ';' '\n' |cut -c1-3)
month_nums=$(seq -w 1 12)

r.proj location=G_ModelE mapset=2x2p5 input=landicefr method=nearest 
r.proj location=G_ModelE mapset=2x2p5 input=landfr method=nearest
r.proj location=G_ModelE mapset=2x2p5 input=ocnfr method=nearest

rasters=$(r.proj location=G_ModelE mapset=2x2p5 -l |grep -E '^temp_|^salt_')
parallel --progress --bar "r.proj -n --q location=G_ModelE mapset=2x2p5 input={1} output={1}_E method=nearest" ::: ${rasters}
parallel --progress --bar "r.grow.distance --q input={1}_E value={1}" ::: ${rasters}

g.region t=${nlevel} -p3
parallel "r.to.rast3 --q input=\$(g.list type=raster pattern=\"{1}_m{2}_z[0-9]?\" sep=,) output={1}_m{2}" ::: temp salt ::: ${month_nums}

# parallel "r3.out.netcdf -p input={1}_m{2} output={1}_m{2}_AQ.nc" ::: temp salt ::: ${month_nums}

# for var in salt temp; do
#   ncrcat -A ${var}_m{01..12}_AQ.nc ${var}_AQ.nc
#   ...
# done

exit

Export to NetCDF

%cd '/home/kdm/projects/GISS/docs' 
import numpy as np
import xarray as xr
import calendar
import glob
import datetime

from grass_session import Session
from grass.script import core as gcore
import grass.script as gscript
# import grass.script.setup as gsetup
# import grass python libraries
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import vector as v
from grass.pygrass.modules.shortcuts import temporal as t
from grass.script import array as garray

afile = glob.glob('JAN*oijl*.nc')[0]
levels = xr.open_dataset(afile)['zoc'].values
months = [_.upper() for _ in list(calendar.month_abbr)[1:]]

for loc in ['GL','AQ']:
    ds = xr.Dataset()
    # S = Session()
    # S.open(gisdb=".", location="G_"+loc, mapset="PERMANENT", create_opts=None)
    with Session(gisdb=".", location="G_"+loc, mapset="PERMANENT", create_opts=None):

        x = garray.array("x", null=np.nan)
        y = garray.array("y", null=np.nan)
        ds['z'] = levels
        ds['x'] = x[0,:]
        ds['y'] = y[:,0]
        ystr = afile[3:7]
        ds['time'] = [datetime.datetime(year=int(ystr), month=m, day=1) for m in range(1,13)]

        for var in ['temp','salt']:
            for i,mon in enumerate(months):
                ii = str(i+1).zfill(2)
                tmp = garray.array3d(f"{var}_m{ii}", null=np.nan)
                ds[f"{var}_m{ii}"] = (('z','y','x'), tmp)

        month_num = [str(m+1).zfill(2) for m in range(len(months))]
        temp_time = xr.concat([ds[f"temp_m{n}"] for n in month_num], dim='time')
        salt_time = xr.concat([ds[f"salt_m{n}"] for n in month_num], dim='time')
        ds['temp'] = (('time','z','y','x'), temp_time.data)
        ds['salt'] = (('time','z','y','x'), salt_time.data)
        for m in month_num:
            ds = ds.drop([f"temp_m{m}",f"salt_m{m}"])

        ds.to_netcdf(f"E2P_{loc}.nc")

/home/kdm/projects/GISS/docs

Display: Antarctic salt at all depths at time 0

ds = xr.open_dataset("E2P_AQ.nc")
_ = ds['salt'].isel({'time':0}).plot(col='z', col_wrap=4, robust=True)

Display: Greenland temp at all times at depth 0

ds = xr.open_dataset("E2P_GL.nc")
_ = ds['temp'].isel({'z':0}).plot(col='time', col_wrap=4, robust=True)

Merge ModelE output with LIME metadata

import numpy as np
import xarray as xr

aq_meta = xr.open_dataset("LIME_AQ.nc")
aq_me = xr.open_dataset("E2P_AQ.nc")
aq = xr.merge([aq_meta,aq_me])

gl_meta = xr.open_dataset("LIME_GL.nc")
gl_me = xr.open_dataset("E2P_GL.nc")
gl = xr.merge([gl_meta,gl_me])

Select appropriate T & S depth for each cell

# # 800 m
# aq = aq.sel(z=800, method="nearest")
# gl = gl.sel(z=800, method="nearest")

# average of top 800 m
aq = aq.sel({'z':slice(0,800)}).mean(dim='z')
gl = gl.sel({'z':slice(0,800)}).mean(dim='z')

NOTDONE Select appropriate T & S for each basin

# aq.groupby('shoreface_basins_100').mean()['salt'].values

Create basin IDs

basins = aq['shoreface_basins_100'].values
for i,b in enumerate(np.unique(basins)):
    if np.isnan(b): basins[np.isnan(basins)] = 0
    basins[basins == b] = i+1
aq['basins'] = (('y','x'), basins)    

basins = gl['shoreface_basins_100'].values
for i,b in enumerate(np.unique(basins)):
    if np.isnan(b): basins[np.isnan(basins)] = 0
    basins[basins == b] = i+1
gl['basins'] = (('y','x'), basins)

Reformat to meet PISM requirements

aq = aq.rename({'temp':'theta_ocean',
                'salt':'salinity_ocean'})
aq['theta_ocean'].attrs['units'] = 'Celsius'
aq['salinity_ocean'].attrs['units'] = 'g/kg'
aq.to_netcdf('ocean_PICO_AQ.nc')

gl = gl.rename({'temp':'theta_ocean',
                'salt':'salinity_ocean'})
gl['theta_ocean'].attrs['units'] = 'Celsius'
gl['salinity_ocean'].attrs['units'] = 'g/kg'
gl.to_netcdf('ocean_PICO_GL.nc')

Test run in PISM

cd ${LIME_ROOT}/runs/ocean.PICO
cp ${LIME_ROOT}/../docs/ocean_PICO_* .

mpiexec -n 4 \
	pismr \
	-i ../pism_Greenland_5km_v1.1.nc \
	-Mx 76 -My 141 -Mz 101 -Mbz 11 \
	-z_spacing equal \
	-Lz 4000 -Lbz 2000 \
	-skip -skip_max 10 \
	-grid.recompute_longitude_and_latitude false \
	-grid.registration corner \
	-surface given \
	-surface_given_file ../pism_Greenland_5km_v1.1.nc \
	-front_retreat_file ../pism_Greenland_5km_v1.1.nc \
	-sia_e 3.0 \
	-stress_balance ssa+sia \
	-topg_to_phi 15.0,40.0,-300.0,700.0 \
	-pseudo_plastic \
	-pseudo_plastic_q 0.5 \
	-till_effective_fraction_overburden 0.02 \
	-tauc_slippery_grounding_lines \
	-ts_file ts_g20km_10ka_hy.nc \
	-extra_file ex_g20km_10ka_hy.nc \
	-extra_vars diffusivity,temppabase,tempicethk_basal,bmelt,tillwat,velsurf_mag,mask,thk,topg,usurf,hardav,velbase_mag,tauc,tendency_of_ice_mass_due_to_discharge,basal_melt_rate_grounded,bmelt \
	-bootstrap \
	-ys 0 -ye 100 \
	-extra_times 0:10:1000 \
	-ts_times 0:yearly:1000 \
	-ocean pico\
	-ocean.pico.file ocean_PICO_GL.nc \
	-ocean.pico.periodic no \
	-frontal_melt.constant.melt_rate 100 \
	-o g20km_1ka_ocean_PICO.nc