notebook.community

Edit and run

James Orr - 3 November 2014

<img align="left" width="40%" src="http://www.lsce.ipsl.fr/Css/img/banniere_LSCE_75.png" >

LSCE/IPSL, CEA-CNRS-UVSQ, Gif-sur-Yvette, France

ORCA025-PISCES model-data comparison: surface zonal means

Synopsis: Makes surface zonal means of DIC and ALK from ORCA025-PISCES and compares them to climatological data.

Procedure: (1) Install and load ferretmagic, (2) read in model output & data files, (3) produce needed ferret "go" scripts (choosebasin, zmsurf, mod2datgrid, dat2refgrid), and (4) use them to plot surface zonal means (model results vs. climatological data).

Plots (see below): zonal mean simulated and observed (1) C$_T$, (2) A$_T$, and (3) A$_T$ $-$ C$_T$.

Install & load ferretmagic extension





In [1]:



    


%install_ext https://raw.github.com/PBrockmann/ipython-ferretmagic/master/ferretmagic.py
%load_ext ferretmagic



    


















    






Installed ferretmagic.py. To use it, type:
  %load_ext ferretmagic

Get MODEL OUTPUT from prodn server on IDRIS http://prodn.idris.fr/thredds/dodsC/IDRISPUBFS/rfry451/ORCA025-PIS2DIC/catalog.html





In [2]:



    


%%ferret -q
  can data/all
  use "http://prodn.idris.fr/thredds/dodsC/IDRISPUBFS/rfry451/ORCA025-PIS2DIC/GRID/mask.nc"
! Problem reading in the following file:
! use "http://prodn.idris.fr/thredds/dodsC/IDRISPUBFS/rfry451/ORCA025-PIS2DIC/MBG/Output/YE/ORCA025-PIS2DIC_20000101_20001231_1Y_ptrc_T.nc"
! use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/IGCM_OUT/ORCA025-PIS2DIC_20000101_20001231_1Y_ptrc_T.nc"
  use "http://prodn.idris.fr/thredds/dodsC/IDRISPUBFS/rfry938/ORCA025-PIS2DIC/MBG/Analyse/TS_MO/ORCA025_1991_2010_1M_DIC2_surf.nc"
  use "http://prodn.idris.fr/thredds/dodsC/IDRISPUBFS/rfry938/ORCA025-PIS2DIC/MBG/Analyse/TS_MO/ORCA025_1991_2010_1M_Alkalini2_surf.nc"

Get DATA from DODS server on ciclad, i.e. http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/gridded





In [3]:



    


%%ferret -q
! GLODAP DIC
  use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/gridded/glodap/TCO2.nc"
! GLODAP ALK
  use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/gridded/glodap/TALK.nc"

Define MYGRID to remove time axis of tmask which is unnecessary and interferes with time_counter axis of model output





In [4]:



    


%%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

choosebasin.jnl - choose basin (mask out all others) using mask file on standard 1°x1° grid)





In [5]:



    


%%writefile choosebasin.jnl
\cancel mode verify

  def sym obasn     $1%|atlarc|ind|pac|indopac|glo<You must specify one of this set: atlarc, ind, pac, indopac, OR glo)%
    
! Read the data (only the first time this routine is called)
  if ($readbmask%|yes>FALSE|*>TRUE%) then
!    Compute Basin mask information (convert from grid start at 20°E to start at 0°)
!    -------------------------------------------------------------------------------
     use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/GRID/LEV_grid.nc"
     use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/specials/LEV_Basins.nc"
     use "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/phase5/DATA/GRID/WOA2001_grid.nc"
!
!    Impose proper axes (instead of just indices in the LEV_Basins file)
     let basin_on_lev20_grid = reshape(basin[d=LEV_Basins.nc],area[d=LEV_grid.nc])
!
!    Put on World Ocean Atlas grid (starting at 0° of longitude)
     let basin_on_lev0_grid = basin_on_lev20_grid[gx=longitude@xact,gy=latitude@xact]
     let bg0 = basin_on_lev0_grid
!
!    Do not consider marginal seas (areas where basin=3)
     let bg  = if bg0 ne 3 then bg0
    
!    Signal that this if BLOCK has been executed (no need to repeat later in same execution of same script)
     def sym readbmask yes
  endif

! Based on user input for obasn (glo, atlarc, ind, pac, indopac), set maskbasin = 1 in that zone (using Levitus mask)
  if ($obasn%|glo>TRUE|*>FALSE%) then
    def sym basntitle "Global"
    def sym basncode "glo"
    let maskbasin = if bg ge -1 then 1
  elif ($obasn%|atlarc>TRUE|*>FALSE%) then
    def sym basntitle "Atlantic"
    def sym basncode "atl"
    let maskbasin = if bg eq 0 or bg0 eq 2 then 1
  elif ($obasn%|pac>TRUE|*>FALSE%) then
    def sym basntitle "Pacific"
    def sym basncode "pac"
    let maskbasin = if bg eq 1 then 1
  elif ($obasn%|ind>TRUE|*>FALSE%) then
    def sym basntitle "Indian"
    def sym basncode "ind"
    let maskbasin = if bg eq 4 then 1
  elif ($obasn%|indopac>TRUE|*>FALSE%) then
    def sym basntitle "Indian + Pacific"
    def sym basncode "indpac"
    let maskbasin = if (bg eq 4) OR (bg eq 1) then 1
  endif

  let maskbasin3D = reshape(maskbasin,mask[d=WOA2001_grid.nc])
        
set mode/last verify



    


















    






Overwriting choosebasin.jnl

zmsurf.jnl - plots surface zonal mean of given variable (must be on rectilinear grid)





In [6]:



    


%%writefile zmsurf.jnl
\cancel mode verify

! Description: plots a zonal mean at surface [k=1] for selected basin
! James Orr, LSCE/IPSL, CEA-CNRS-UVSQ, CEA Saclay, 1 Nov 2014
! ============================================================

  def sym var     $1%CO3%
  def sym basn    $2%glo%
  def sym xlabl   $3%no%
  def sym ylabl   $4%no%
  def sym latmin  $5%+70%
  def sym latmax  $6%+90%
  def sym pltopts $7
  def sym overplt $8%over%

! Select basin (mask out all others)
  go choosebasin ($basn)
    
  let ($var)bmskd  = ($var) * maskbasin3D

  if ($overplt%|over>FALSE|*>TRUE%) then  !Pre-plot setup (only if it is the 1st plot, i.e., not an overplot)
!   ppl axnmtc 4, 4 
!   ppl axlint 1, 1
    ppl labset 0.2,0.19,0.19,0.15
    ppl axlsze 0.19, 0.19
!   ppl axlabp +1, -1
  endif

  plot/($pltopts) ($var)bmskd[x=@ave,k=1]

  if ($overplt%|over>FALSE|*>TRUE%) then  !Post-plot setup (only if it is the 1st plot; not an overplot)
!   ppl xfor ($xformat)
    ppl yfor ($yformat)
!   ppl cross 1
    if ($xlabl%|no>FALSE|*>TRUE%) then
      ppl xlab ($xlabl)
    endif
    if ($ylabl%|no>FALSE|*>TRUE%) then
      ppl ylab ($ylabl)
    endif
    ppl plot
  endif

set mode/last verify



    


















    






Overwriting zmsurf.jnl

mod2datgrid.jnl - routine to interpolate model output to data grid (regular grid)





In [7]:



    


%%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) griddat
! sho grid griddat
    
  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



    


















    






Overwriting mod2datgrid.jnl

dat2refgrid.jnl - routine to interpolate data to reference grid (both are rectilinear grids)





In [8]:



    


%%writefile dat2refgrid.jnl
\cancel mode verify

  def sym vardat   $1%ALK%
  def sym varref   $2%mask%
    
  define grid/like=($varref) gridref
    
  set mode interpolate
  let ($vardat)_refgrid = ($vardat)[gx=gridref, gy=gridref] 

set mode/last verify



    


















    






Overwriting dat2refgrid.jnl

Define proper lat, lon coordinates from MASK and MYGRID





In [9]:



    


%%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

Plot 1: C$_T$ zonal means (ORCA025 and data)

Prepare zonal mean DIC model output & data (put on same regular grid)





In [10]:



    


%%ferret -q
! Get basin masks and reference grid (World Ocean Atlas, 1° x 1°)
  def sym readbmask "notyet"
  go choosebasin glo

! DIC model output:
! For tests, use just 1 timestep (for final product (use decadal average 1990-1999
! - namely: DIC[d=2,t=01-jan=1990:31-dec-1999@ave]
! - better still would be to read in the decadal average directly (different input file)
  let var=DIC2[d=2,t=1-jul-1994]/1.028e-6
! let myvar=if TMASK[d=mask.nc,g=mygrid] eq 1 then var   !This does not quite work (come back to this later)
  let myvar = if var then var

! Specify model field (datm) and data field (datd)
  let datm = myvar[k=1]
  let datd = TCO2[d=TCO2.nc,k=1]

! Specify reference field (the grid onto which both the model & the data will be placed)
  let ref = mask[d=WOA2001_grid.nc, k=1]

! Regrid model output and data to reference grid (World Ocean Atlas)
  go mod2datgrid datm ref lon_masked lat_masked
  let datm_ref = datm_reg

! Regrid data to reference grid (World Ocean Atlas) - produces "datd_refgrid"
  go dat2refgrid datd ref
  let datd_ref = reshape(datd_refgrid, datm_ref)

Plot zonal means (Global, Atlantic, Indian, Pacific)





In [11]:



    


%%ferret -q -s 800,400
!%%ferret -q -p -f zmsurf_ORCA025-PISCES_DIC.pdf

! Choose color palette & define line patterns:
  def sym opal nair1_cmyk
! def sym opal light_centered

  def sym dot "(.025,.025,.025,.025)"
  def sym dotwide "(.045,.045,.045,.045)"
  def sym dashdot "(.1,.05,.025,.05)"
  def sym dotwider "(.045,.1,.045,.1)"

! Make 3 plots:
! -------------
  go page_new -l 2 2 5 95 10 90 0 0
  set win/clear
  go boldf_plattner

! Refine plots with a some ferret (PPL) setup commands
  ppl color 3    0, 55,   0    !dark green instead of default green (too light)
  ppl axlabp -1, -1            !labeled axes on left and bottom
  ppl axnmtc 2, 4              !number of minor ticks on x and y axes
  ppl axlint 1, 1              !numer  of major ticks to label (1 = every, 2 = every other)
  def sym yformat (i5)

! Specify y axis limits 
  def sym ymin 1800 ; def sym ymax 2200 ; def sym yinc 100
  def sym ylims ($ymin):($ymax):($yinc)
  let ysubtitle = `($ymin)` + 0.9*`($ymax) - ($ymin)`

! 1) Global
  set view 1
  go margins_set 15 8 13 5
  go zmsurf datm_ref glo no "C_T (($greek)m($roman)mol kg^-^1)"  -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref glo no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Global"

! 2) Atlantic
  set view 2
  go margins_set 15 8 18 0
  go zmsurf datm_ref atlarc no no -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref atlarc no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Atlantic"

! 3) Indian
  set view 3
  go margins_set 15 8 13 5
  go zmsurf datm_ref ind no "C_T (($greek)m($roman)mol kg^-^1)"  -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref ind no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Indian"

! 4) Pacific
  set view 4
  go margins_set 15 8 18 0
  go zmsurf datm_ref pac no no -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref pac no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Pacific"

Fig. 1. Surface zonal means of C$_T$ ($\mu$mol kg$^{-1}$) for the GLODAP data (blue) and ORCA025-PISCES (red) over the global ocean as well as the Atlantic, Indian, and Pacific basins.

ORCA025-PISCES generally underestimates the GLODAP zonal mean C$_T$ particularly in the lower latitudes of the the Atlantic and Indian Oceans. Global zonal means for the model are low by up to 40 $\mu$mol kg$^{-1}$.

Plot 2: A$_T$ zonal means (ORCA025 and data)

Prepare zonal mean ALK model output & data (put on same regular grid)





In [12]:



    


%%ferret -q
! Alk model output:
  let var=Alkalini2[d=3,t=1-jul-1994]/1.028e-6
! let myvar=if TMASK[d=mask.nc,g=mygrid] eq 1 then var   !This does not quite work (come back to this later)
  let myvar = if var then var

! Specify model field (datm) and data field (datd)
  let datm = myvar[k=1]
  let datd = TALK[d=TALK.nc,k=1]
! Regrid model output and data to reference grid (World Ocean Atlas)
  go mod2datgrid datm ref lon_masked lat_masked
  let datm_ref = datm_reg

! Regrid data to reference grid (World Ocean Atlas) - produces "datd_refgrid"
  go dat2refgrid datd ref
  let datd_ref = reshape(datd_refgrid, datm_ref)

Plot zonal means (Global, Atlantic, Indian, Pacific)





In [13]:



    


%%ferret -q -s 800,400
!%%ferret -q -p -f zmsurf_ORCA025-PISCES_DIC.pdf
! Choose color palette & define line patterns:
  def sym opal nair1_cmyk
! def sym opal light_centered
! Make 3 plots:
! -------------
  go page_new -l 2 2 5 95 10 90 0 0
  set win/clear
  go boldf_plattner
! Refine plots with a some ferret (PPL) setup commands
  ppl axlabp -1, -1            !labeled axes on left and bottom
  ppl axnmtc 2, 4              !number of minor ticks on x and y axes
  ppl axlint 1, 1              !numer  of major ticks to label (1 = every, 2 = every other)
  def sym yformat (i5)

! Specify y axis limits 
  def sym ymin 2050 ; def sym ymax 2450 ; def sym yinc 100
  def sym ylims ($ymin):($ymax):($yinc)
  let ysubtitle = `($ymin)` + 0.9*`($ymax) - ($ymin)`

! 1) Global
  set view 1
  go margins_set 15 8 13 5
  go zmsurf datm_ref glo no "A_T (($greek)m($roman)mol kg^-^1)"  -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref glo no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Global"

! 2) Atlantic
  set view 2
  go margins_set 15 8 18 0
  go zmsurf datm_ref atlarc no no -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref atlarc no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Atlantic"

! 3) Indian
  set view 3
  go margins_set 15 8 13 5
  go zmsurf datm_ref ind no "A_T (($greek)m($roman)mol kg^-^1)"  -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref ind no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Indian"

! 4) Pacific
  set view 4
  go margins_set 15 8 18 0
  go zmsurf datm_ref pac no no -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref pac no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Pacific"

Fig. 2. Surface zonal means of A$_T$ ($\mu$mol kg$^{-1}$) for the GLODAP data (blue) and ORCA025-PISCES (red) over the global ocean as well as the Atlantic, Indian, and Pacific basins.

Patterns of differences between ORCA025-PISCES and the GLODAP data for A$_T$ are generally consistent with those for C$_T$. Once again the model underestimates observed values almost everywhere. The largest discrepencies are found in the Atlantic basin and in the low latitudes of the Indian basin.

Plot 3: A$_T$ $-$ C$_T$ zonal means (ORCA025 and data)

Prepare zonal mean ALK-DIC model output & data (put model difference & data difference on same regular grid)





In [14]:



    


%%ferret -q
! Alk - DIC model output:
  let var = (Alkalini2[d=3,t=1-jul-1994] - DIC2[d=2,t=1-jul-1994])/1.028e-6
! let myvar=if TMASK[d=mask.nc,g=mygrid] eq 1 then var   !This does not quite work (come back to this later)
  let myvar = if var then var

! Specify model field (datm) and data field (datd)
  let datm = myvar[k=1]
  let datd = TALK[d=TALK.nc,k=1] - TCO2[d=TCO2.nc,k=1]
! Regrid model output and data to reference grid (World Ocean Atlas)
  go mod2datgrid datm ref lon_masked lat_masked
  let datm_ref = datm_reg

! Regrid data to reference grid (World Ocean Atlas) - produces "datd_refgrid"
  go dat2refgrid datd ref
  let datd_ref = reshape(datd_refgrid, datm_ref)

Plot zonal means (Global, Atlantic, Indian, Pacific)





In [15]:



    


%%ferret -q -s 800,400
!%%ferret -q -p -f zmsurf_ORCA025-PISCES_DIC.pdf
! Choose color palette & define line patterns:
  def sym opal nair1_cmyk
! def sym opal light_centered
! Make 3 plots:
! -------------
  go page_new -l 2 2 5 95 10 90 0 0
  set win/clear
  go boldf_plattner
! Refine plots with a some ferret (PPL) setup commands
  ppl axlabp -1, -1            !labeled axes on left and bottom
  ppl axnmtc 2, 4              !number of minor ticks on x and y axes
  ppl axlint 1, 1              !numer  of major ticks to label (1 = every, 2 = every other)
  def sym yformat (i5)
! Specify y axis limits 
  def sym ymin 100 ; def sym ymax 400 ; def sym yinc 100
  def sym ylims ($ymin):($ymax):($yinc)
  let ysubtitle = `($ymin)` + 0.9*`($ymax) - ($ymin)`

! 1) Global
  set view 1
  go margins_set 15 8 13 5
  go zmsurf datm_ref glo no "A_T - C_T (($greek)m($roman)mol kg^-^1)"  -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref glo no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Global"

! 2) Atlantic
  set view 2
  go margins_set 15 8 18 0
  go zmsurf datm_ref atlarc no no -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref atlarc no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Atlantic"

! 3) Indian
  set view 3
  go margins_set 15 8 13 5
  go zmsurf datm_ref ind no "A_T - C_T (($greek)m($roman)mol kg^-^1)"  -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref ind no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Indian"

! 4) Pacific
  set view 4
  go margins_set 15 8 18 0
  go zmsurf datm_ref pac no no -90 +90 "/nolab/set/col=red/dash/thick=2/hlim=90s:90n:30/vlim=($ylims)/line" firstplt
  go zmsurf datd_ref pac no no -90 +90 "/nolab/over/col=blue/thick=2/dash=($dotwide)  " over
  label 87 `ysubtitle` +1 0 .20 "Pacific"

Fig. 3. Surface zonal means of A$_T$ $-$ C$_T$ ($\mu$mol kg$^{-1}$) for the GLODAP data (blue) and ORCA025-PISCES (red) over the global ocean as well as the Atlantic, Indian, and Pacific basins.

Differences between ORCA025-PISCES and the GLODAP data for A$_T$ $-$ C$_T$ are much smaller than those for either A$_T$ and C$_T$. That is surface C$_T$ tends to equlibrate with the atmosphere via air-sea CO$_2$ exchange and the chemical capacity of the surface ocean to absorbe CO$_2$ is determined by the A$_T$ $-$ C$_T$ difference.

Content source: jamesorr/Atlas_orca025-pisces

Similar notebooks:

notebook.community | gallery | about