<img align="left" width="40%" src="http://www.lsce.ipsl.fr/Css/img/banniere_LSCE_75.png" >
Make surface maps of some key biogeochemical variables from ORCA025-PISCES and compare to climatological data. To do this very basic, qualitative evaluation, use pyferret in IPython Notebook.
After installing ferretmagic and loading it, this Notebook (1) makes a list of files to be treated, (2) makes and writes out three new plotting routines (Ferret go scripts: mapmod.jnl, mapdat.jnl, and mod2datgrid.jnl), and (3) uses those routines in subsequent cells to make surface maps of model results vs. climatological data.
In [3]:
%install_ext https://raw.github.com/PBrockmann/ipython-ferretmagic/master/ferretmagic.py
In [4]:
%load_ext ferretmagic
In [5]:
%%ferret -q
can data/all
use "http://prodn.idris.fr/thredds/dodsC/IDRISPUBFS/rfry451/ORCA025-PIS2DIC/GRID/mask.nc"
use "http://prodn.idris.fr/thredds/dodsC/IDRISPUBFS/rfry451/ORCA025-PIS2DIC/MBG/Output/YE/ORCA025-PIS2DIC_20000101_20001231_1Y_ptrc_T.nc"
In [6]:
%%ferret -q
! CARS phosphate
! use "/home/biomac1/blevu/DATA/CARS2009/phosphate_cars2009.nc"
use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/gridded/CARS2009/phosphate_cars2009.nc"
! use "/home/biomac2/jpalm/WOA/phosphate_annual_1deg.nc"
! CARS nitrate
! use "/home/biomac1/blevu/DATA/CARS2009/nitrate_cars2009.nc"
use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/gridded/CARS2009/nitrate_cars2009.nc"
! use "/home/biomac2/jpalm/WOA/nitrate_annual_1deg.nc"
! GLODAP DIC
! use "/home/biomac1/blevu/DATA/GLODAP/tco2_ann_Key.cdf"
use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/gridded/glodap/TCO2.nc"
! GLODAP ALK
! use "/home/biomac1/blevu/DATA/GLODAP/talk_ann_Key.cdf"
use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/gridded/glodap/TALK.nc"
In [7]:
%%ferret -q
def grid/x=x/y=y/z=z mygrid
set axis/name=xlon X
set axis/name=ylat Y
set axis/name=zdepth Z
In [8]:
%%writefile mapmod.jnl
\cancel mode verify
def sym var $1%DIC%
def sym lon $2%lon%
def sym lat $3%lat%
def sym xlabl $4%no%
def sym ylabl $5%no%
def sym pltopts $6
! James Orr, LSCE/IPSL, CEA-CNRS-UVSQ, CEA Saclay, 15 Oct 2014
! ============================================================
! let varl = ($var)[x=@ave,y=($latmin):($latmax)@ave]
let varl = ($var)
ppl labset 0.2,0.16,0.16,0.15
ppl axlsze 0.17, 0.17
ppl axnmtc 2, 2
ppl axlint 2, 1
ppl axlabp -1, -1
shade/nolab/set/hlim=120e:480e:30/vlim=90s:90n:30/mod/($pltopts) ($var), ($lon), ($lat)
if ($xlabl%|no>FALSE|*>TRUE%) then
ppl xlab ($xlabl)
endif
if ($ylabl%|no>FALSE|*>TRUE%) then
ppl ylab ($ylabl)
endif
ppl shade
go fland 60 gray36; go land
set mode/last verify
In [9]:
%%writefile mapdat.jnl
\cancel mode verify
def sym var $1%DIC%
def sym xlabl $2%no%
def sym ylabl $3%no%
def sym pltopts $4
! James Orr, LSCE/IPSL, CEA-CNRS-UVSQ, CEA Saclay, 15 Oct 2014
! ============================================================
! let varl = ($var)[x=@ave,y=($latmin):($latmax)@ave]
let varl = ($var)
ppl labset 0.2,0.16,0.16,0.15
ppl axlsze 0.17, 0.17
ppl axnmtc 2, 2
ppl axlint 2, 1
ppl axlabp -1, -1
! Define as a module axis (e.g., not set in CARS) needed to show full map from 120E to 120E
set axis/mod `($var),RETURN=xaxis`
shade/nolab/set/x=240w:120e/hlim=120:480:30/vlim=90s:90n:30/mod/($pltopts) ($var)
if ($xlabl%|no>FALSE|*>TRUE%) then
ppl xlab ($xlabl)
endif
if ($ylabl%|no>FALSE|*>TRUE%) then
ppl ylab ($ylabl)
endif
ppl shade
go fland 60 gray36; go land
set mode/last verify
In [10]:
%%writefile mod2datgrid.jnl
\cancel mode verify
def sym varmod $1%DIC%
def sym vardat $2%ALK%
def sym lonsmod $3
def sym latsmod $4
define grid/like=($vardat) gridref
! sho grid gridref
let lonsd = x[gx=`($vardat),return=xaxis`]
let latsd = y[gy=`($vardat),return=yaxis`]
let latsm = xsequence(($latsmod))
let lonsm = xsequence(($lonsmod))
! Fix problem when model longitude (range -180 to +180) is outside data longitude range (0 to 360.)
let lonsdmin = lonsd[i=@min]
let lonsdmax = lonsd[i=@max]
let lonsm_lo = if lonsm lt lonsdmin then lonsm + 360. else lonsm
let lonsmd = if lonsm_lo gt lonsdmax then lonsm_lo - 360. else lonsm_lo
let ($varmod)_reg = scat2grid_bin_xy(lonsmd, latsm, xsequence(($varmod)), lonsd, latsd)
set mode/last verify
In [11]:
%%ferret
let lon_masked=if TMASK[g=mygrid, d=mask.nc, k=1] eq 1 then nav_lon[d=mask.nc]
let lat_masked=if TMASK[g=mygrid, d=mask.nc, k=1] eq 1 then nav_lat[d=mask.nc]
set variable/units="degrees_east"/title="Longitude" lon_masked
set variable/units="degrees_north"/title="Latitude" lat_masked
In [12]:
%%ferret
! Consider only every 4th data point in x and y directions to conserve computational time (roughly 1/4° -> 1°)
def sym subreg=i=1:1442:4,j=1:1021:4
! Consider every data point in x and y directions (more costly; finer details may not be noticible on globap maps)
! def sym subreg=i=1:1442:1,j=1:1021:1
! set mem/siz=300
To use the custom palettes below, you may need to add my directory to your FER_PALETTE environment variable. If the output below does not end in "/home/users/orr/ferret/palettes", you can manually do this command in your .bashrc file or make the equivalent command in your .csrhc file (when running csh or tcsh).
In [13]:
%%bash
echo $FER_PALETTE
After adding the above directory to the FER_PALETTE environment variable, you will have access to the color palettes used below.
In [14]:
%%ferret -q -s 800,400
!%%ferret -q -p -f surfmap_ORCA025-PISCES_PO4.pdf
set mode verify
let var=PO4[d=2]/122*1e6
let myvar=if TMASK[d=mask.nc,g=mygrid] eq 1 then var
! Choose color palette: most below are custom made (put asterix:/home/users/orr/ferret/palettes in $FER_PALETTE env variable)
! def sym opal rain_cmyk
! def sym opal rain0_cmyk
! def sym opal lodyc
! def sym opal myblue_darkorange
def sym opal nair1_cmyk
! def sym opal light_centered
! specify model field (datm) and data field (datd)
let datm = myvar[k=1,l=1]
let datd = mean[d=phosphate_cars2009.nc,k=1]
! Make 3 plots:
! -------------
go page_new 1 3 5 95 0 100 0 0
set win/clear
go boldf_plattner
! 1) Model
set view 1
go margins_set 15 8 8 10
go mapmod datm[($subreg)] lon_masked[($subreg)] lat_masked[($subreg)] no Latitude "/lev=(0,2,0.2)/pal=($opal)"
! 2) Data
set view 2
go margins_set 8 15 8 10
go mapdat datd no Latitude "/lev=(0,2,0.2)/pal=($opal)"
! 3) Model - Data difference
set view 3
go margins_set 1 22 8 10
def sym opal light_centered
go mod2datgrid datm datd lon_masked lat_masked
! Compute model - data difference on data grid
! -> reshape command below works around bug in Ferret 6.9 (should not be needed; without it, cannot make map w/ right LON limits)
let diffmd = reshape(datm_reg - datd, datd)
go mapdat diffmd Longitude Latitude "/lev=(-inf)(-1.4,1.4,0.2)(-0.1,0.1,0.2)(inf)/pal=($opal)"
ORCA025-PISCES overestimates observed surface phosphate concentrations in the low latitudes and underestimates the same in the high latitudes, particularly in HNLC regions (North Pacific and Southern Ocean). In the Arctic, the model has particularly high values in coastal waters just north of Siberia; conversely, available observations (summer biased) indicate particularly low values. The largest differences appear to be in or near ice covered regions. Model results in marginal seas may stand out visually but are not to be taken seriously due in part to poor initial conditions as well as inadequate resolution and tuning.
In [15]:
%%ferret -q -s 800,400
!%%ferret -q -p -f surfmap_ORCA025-PISCES_NO3.pdf
let var=NO3[d=2]*16/122*1e6
let myvar=if TMASK[d=mask.nc,g=mygrid] eq 1 then var
let datm = myvar[k=1,l=1]
let datd = mean[d=nitrate_cars2009.nc,k=1]
! let datd = n_an[d=nitrate_annual_1deg.nc,k=2]
go page_new 1 3 5 95 0 100 0 0
set win/clear
go boldf_plattner
! 1) Model
set view 1
go margins_set 15 8 8 10
def sym opal nair1_cmyk
go mapmod datm[($subreg)] lon_masked[($subreg)] lat_masked[($subreg)] no Latitude "/lev=(0,30,3)(O.5,1,0.5)/pal=($opal)"
! 2) Data
set view 2
go margins_set 8 15 8 10
go mapdat datd no Latitude "/lev=(0,30,3)(0.5,1,0.5)/pal=($opal)"
! 3) Model - Data difference
set view 3
go margins_set 1 22 8 10
def sym opal light_centered
! Regrid model to data grid
go mod2datgrid datm datd lon_masked lat_masked
! Compute model - data difference (on data grid)
! -> reshape command below works around bug in Ferret 6.9 (should not be needed; without it, cannot make map w/ right LON limits)
let diffmd = reshape(datm_reg - datd, datd)
! Make map
go mapdat diffmd Longitude Latitude "/lev=(-inf)(-24,24,4)(-2,2,4)(inf)/pal=($opal)"
As seen in Fig. 1 for phosphate, ORCA025-PISCES also underestimates surface nitrate in the high latitudes, particularly in HNLC regions (North Pacific and Southern Ocean). In the low latitudes model-data differences are generally small in absolute terms. Larger absolute differences are seen in regions of mode water formation in the Southern Ocean and North Pacific. In the Arctic, simulated nitrate appears much larger than observed in coastal waters just north of Siberia, as seen for phosphate (Fig. 1). Elsewhere in the Arctic simulated surface nitrate is generally larger than observed, whereas simulated surface phosphate is generally less than observed.
In [16]:
%%ferret -q -s 800,400
!%%ferret -q -p -f surfmap_ORCA025-PISCES_Ct.pdf
let var=DIC2[d=2]/1.028e-6
let myvar=if TMASK[d=mask.nc,g=mygrid] eq 1 then var
let datm = myvar[k=1,l=1]
let datd = TCO2[d=TCO2.nc,k=2]
go page_new 1 3 5 95 0 100 0 0
set win/clear
go boldf_plattner
! 1) Model
set view 1
go margins_set 15 8 8 10
def sym opal nair1_cmyk
go mapmod datm[($subreg)] lon_masked[($subreg)] lat_masked[($subreg)] no Latitude "/lev=(-inf)(1800,2300,25)(inf)/pal=($opal)"
! 2) Data
set view 2
go margins_set 8 15 8 10
go mapdat datd no Latitude "/lev=(-inf)(1800,2300,25)(inf)/pal=($opal)"
! 3) Model - Data difference
set view 3
go margins_set 1 22 8 10
def sym opal light_centered
! Regrid model to data grid
go mod2datgrid datm datd lon_masked lat_masked
! Compute model - data difference (on data grid)
! -> reshape command below works around bug in Ferret 6.9 (should not be needed; without it, cannot make map w/ right LON limits)
let diffmd = reshape(datm_reg - datd, datd)
! Make map
go mapdat diffmd Longitude Latitude "/lev=(-inf)(-100,100,20)(-10,10,5)(inf)/pal=($opal)"
Simulated surface C$_T$ in ORCA025-PISCES overestimates the observed GLODAP climatology in the subarctic Pacific and parts of the equatorial Pacific and Southern Ocean (by up to 60 $\mu$mol kg$^{-1}$), while it underestimates the GLODAP reference by at least as much in the Northern Indian Ocean as well as the low and mid latitudes of the Atlantic Ocean. The pattern of high C$_T$ in coastal waters just north of Siberia is similar to that for nutrient concentrations (Figs. 1 and 2).
In [17]:
%%ferret -q -s 800,400
!%%ferret -q -p -f surfmap_ORCA025-PISCES_At.pdf
let var=ALKALINI2[d=2]/1.028e-6
let myvar=if TMASK[d=mask.nc,g=mygrid] eq 1 then var
let datm = myvar[k=1,l=1]
let datd = TALK[d=TALK.nc,k=1]
go page_new 1 3 5 95 0 100 0 0
set win/clear
go boldf_plattner
! 1) Model
set view 1
go margins_set 15 8 8 10
def sym opal nair1_cmyk
go mapmod datm[($subreg)] lon_masked[($subreg)] lat_masked[($subreg)] no Latitude "/lev=(-inf)(1900,2400,25)(inf)/pal=($opal)"
! 2) Data
set view 2
go margins_set 8 15 8 10
go mapdat datd no Latitude "/lev=(-inf)(1900,2400,25)(inf)/pal=($opal)"
! 3) Model - Data difference
set view 3
go margins_set 1 22 8 10
def sym opal light_centered
! Regrid model to data grid
go mod2datgrid datm datd lon_masked lat_masked
! Compute model - data difference (on data grid)
! -> reshape command below works around bug in Ferret 6.9 (should not be needed; without it, cannot make map w/ right LON limits)
let diffmd = reshape(datm_reg - datd, datd)
! Make map
go mapdat diffmd Longitude Latitude "/lev=(-inf)(-100,100,20)(-10,10,5)(inf)/pal=($opal)"
The difference between simulated surface A$_T$ and the observed GLODAP climatology resembles that for C$_T$ (Fig. 3), namely with underestimates by more than 80 $\mu$mol kg$^{-1}$ in the Northern Indian Ocean as well as the low and mid latitudes of the Atlantic Ocean. Yet overestimates of simulated A$_T$ are smaller in the equatorial Pacific and Southern Ocean (less than 20 $\mu$mol kg$^{−1}$). In Siberian coastal waters, we see the same pattern of high A$_T$ as seen for C$_T$ and nutrients (Figs. 1-3).
In [18]:
%%ferret -q -s 800,400
!%%ferret -q -p -f surfmap_ORCA025-PISCES_At-Ct.pdf
set mem/siz=300
let var=(ALKALINI2[d=2] - DIC[d=2])/1.028e-6
let myvar=if TMASK[d=mask.nc,g=mygrid] eq 1 then var
let datm = myvar[k=1,l=1]
let datd = TALK[d=TALK.nc,k=2] - TCO2[d=TCO2.nc,k=2]
go page_new 1 3 5 95 0 100 0 0
set win/clear
go boldf_plattner
def sym olevs (-inf)(120,380,20)(inf)
! 1) Model
set view 1
go margins_set 15 8 8 10
def sym opal nair1_cmyk
go mapmod datm[($subreg)] lon_masked[($subreg)] lat_masked[($subreg)] no Latitude "/lev=($olevs)/pal=($opal)"
! 2) Data
set view 2
go margins_set 8 15 8 10
go mapdat datd no Latitude "/lev=($olevs)/pal=($opal)"
! 3) Model - Data difference
set view 3
go margins_set 1 22 8 10
def sym opal light_centered
! Regrid model to data grid
go mod2datgrid datm datd lon_masked lat_masked
! Compute model - data difference (on data grid)
! -> reshape command below works around bug in Ferret 6.9 (should not be needed; without it, cannot make map w/ right LON limits)
let diffmd = reshape(datm_reg - datd, datd)
! Make map
go mapdat diffmd Longitude Latitude "/lev=(-inf)(-50,50,10)(-10,10,5)(inf)/pal=($opal)"
The simulated difference in A$_T$ $-$ C$_T$ may be used for a very rough proxy of carbonate ion concentration (a better proxy would be A$_C$ $-$ C$_T$, but A$_C$ would have to be calculated offline). The model-data discrepancy for A$_T$ $-$ C$_T$ difference is usually much smaller than seen for the A$_T$ and C$_T$ tracers individually. The model overestimates this difference in Eastern Boundary Upwelling Systems (EBUS), e.g., in the California Current, Humboldt Current, and Canary Current systems. It also overestimates the same quantity offshore of Southeast Asia, Eastern Australia, and in the Arabian Sea. However in these near coastal areas, GLODAP data coverage is sparse at best. The largest underestimates of the A$_T$ $-$ C$_T$ difference are found in Sea of Okhotsh and Bering Sea, the subpolar Northwestern Atlantic, and the Bay of Bengal. Over most of the open ocean, the simulated A$_T$ $-$ C$_T$ difference overestimates the GLODAP based estimate by typically about 20 $\mu$mol kg$^{-1}$. The opposite tendancy is also found in the opean ocean, e.g., underestimates by around 20 $\mu$mol kg$^{-1}$) occurs over much of the Southern Ocean, namely over most of the Atlantic and Indian sectors as well as about one-third of the Pacific sector.