source: trunk/src/fig_mld.pro @ 10

;+
;
; ===========
; fig_mld.pro
; ===========
;
; .. program:: fig_mld
;
; main POMME project program
;
; - read the 2 mld models :file:`${PROJECT_ID}/noeddy_mld.nc` and :file:`${PROJECT_ID}/eddy_mld.nc`
; - read mld insitu data and metadata
; - compute daily mean insitu mld data
; - compute daily mean modeled data in boxes of 1/2° and 1° around insitu
;   measurement location
; - plot +todo+
; - write :file:`${PROJECT_OD}/fig_mld.txt` containing for each day between
;   20000927 and +todo+
;
;   * time step (0 is 20000927)
;   * latitude
;   * longitude
;   * insitu daily mean of MLD critere 01
;   * insitu daily mean of MLD critere 02
;   * insitu daily mean of MLD critere 05
;   * insitu daily mean of MLD critere 10
;   * insitu daily mean of MLD critere melch
;   * model 1 mean of MLD in the box latidude+-0.5°,longitude+-.5°
;   * model 2 mean of MLD in the box latidude+-0.5°,longitude+-.5°
;   * model 1 stddev of MLD in the box latidude+-0.5°,longitude+-.5°
;   * model 2 stddev of MLD in the box latidude+-0.5°,longitude+-.5°
;   * model 1 mean of MLD in the box latidude+-1°,longitude+1°
;   * model 2 mean of MLD in the box latidude+-1°,longitude+1°
;   * model 1 stddev of MLD in the box latidude+-1°,longitude+1°
;   * model 2 stddev of MLD in the box latidude+-1°,longitude+1°
;
; .. only:: man
;
;    Figure is visible on PDF and HTML documents only.
;
; .. only:: html
;
;     .. graphviz::
;
;        digraph fig_mld {
;           graph [
;           rankdir="TB",
;           ]
;           ;
;           pomme_asc [shape=ellipse,fontname=Courier,label="${VARAMMA_ID}/*.asc",
;           URL="http://www.locean-ipsl.upmc.fr/~fplod/varamma/varamma_ws/doc/html/sphinx/guides/data_content.html#data_station"];
;           pomme_prn [shape=ellipse,fontname=Courier,label="${VARAMMA_ID}/*.prn",
;           URL="http://www.locean-ipsl.upmc.fr/~fplod/varamma/varamma_ws/doc/html/sphinx/guides/data_content.html#data_mld_insitu"];
;           pomme_nc [shape=ellipse,fontname=Courier,label="${VARAMMA_ID}/*.nc",
;           URL="http://www.locean-ipsl.upmc.fr/~fplod/varamma/varamma_ws/doc/html/sphinx/guides/data_content.html#data_mld_models"];
;           pomme_ps [shape=ellipse,fontname=Courier,label="${VARAMMA_OD}/fig_mld_valid_sta.ps"];
;
;           fig_mld [shape=box,
;           fontname=Courier,
;           color=blue,
;           URL="http://forge.ipsl.jussieu.fr/pomme/browser/trunk/src/fig_mld.pro"
;           label="${PROJECT}/src/fig_mld.pro"];
;
;           {pomme_asc pomme_prn pomme_nc} -> {fig_mld} -> {pomme_ps};
;        }
;
; :examples:
;
; .. code-block:: guess
;
;    IDL> @init_pom
;    IDL> fig_mld
;    IDL> exit
;
; :file:`${PROJECT_OD}/fig_mld_valid_sta.ps` can be used for future work.
;
; :file:`${PROJECT_OD}/fig_mld.txt` can be used for future work (mainly
; comparison).
;
; :uses:
;
; :func:`read_ncdf <saxo:read_ncdf>`
; :func:`julday <saxo:julday>`
; :func:`plt1d <saxo:plt1d>`
; :func:`moyenne <saxo:moyenne>`
; :func:`domdef <saxo:domdef>`
;
; :func:`file_prn_to_mem`
; :func:`file_asc_to_mem`
; :func:`dmean_mld`
; :func:`join_dmean_asc`
;
; :func:`reform`
;
; :see also:
;
; figure 4 of submitted paper +todo+
;
; :ref:`project_init.pro`
; :ref:`cm_project.pro`
;
; :todo:
;
; 
; - fplod 20110426T170242Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * make it work : now pb with ::
;
;     % Error: impossible to make a "t" type plot with jpt = 1
;     % Attempt to subscript MLD1 with IDAY is out of range.
;
; remove hard coded value (95) of time axe of model correction
;
; better organisation of output files
; ascii output file
;
; 250 = red
; 90 = light blue
; 150 = green
; ++ longitude should be dark blue
; ++ latitude should be orange
;
; superposer plot latitude et longitude
;
; mettre les labels et le titre des Y du plot longitude à droite
; see x recipe http://www.dfanning.com/graphics_tips/xattop.html
;
; how to avoid lines between legs without retouche++
;
; :history:
;
; - fplod 20110426T161700Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * rely on project_init.pro instead of init_pom.pro
;   * fix typo (copy/paste std1 vs mld5)
;   * add init_pom
;   * handle 2D vs 1D my* because ml use only delta = 0
;
; - fplod 20101209T100325Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * align ytitles to final figures
;   * add ascii output
;
; - fplod 20101208T140038Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * revision according to Marina/myself work (20101130)
;     tricky convention differences on longitude between in situ and model files
;   * revision according to Marina work (20101201)
;     plots production, selection of delta = 1 for figure
;   * add graphviz
;
; - fplod 20101129T151939Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * loop on delta
;   * create output file
;   * use http://www.dfanning.com/math_tips/nans.html to test *equality* to NaN
;
; - fplod 20101129T130822Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * add reading of insitu MLD data files and associated
;   * use the right filenames for models
;   * complete header
;
; - fplod 20101125T145638Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * add ReST minimal header
;
; - marina 20101125T130000Z
;
;   * creation draft
;
;-
PRO fig_mld
;
compile_opt idl2, strictarrsubs
;
@cm_project
;
; read MLD in insitu files and metafiles
data_asc=file_asc_to_mem()
data_prn=file_prn_to_mem()
data_dmean=dmean_mld(data_prn)
data_insitu=join_dmean_asc(data_dmean,data_asc)
;
; temps des données insitu
time_insitu=data_insitu.doy
;
; build compatible insitu arrays with models arrays
ndays = n_elements(time_insitu)
mld1 =  findgen(ndays)
mld2 = findgen(ndays)
mld3 = findgen(ndays)
mld4 = findgen(ndays)
mld5 = findgen(ndays)
std1 = findgen(ndays)
xlon = findgen(ndays)
ylat = findgen(ndays)
FOR iday=0, ndays - 1 DO BEGIN
   indice_insitu=WHERE(time_insitu-time_insitu[0] EQ iday,count)
   IF (count EQ 0) THEN BEGIN
      mld1[iday] = !VALUES.F_NAN
      mld2[iday] = !VALUES.F_NAN
      mld3[iday] = !VALUES.F_NAN
      mld4[iday] = !VALUES.F_NAN
      mld5[iday] = !VALUES.F_NAN
      std1[iday] = !VALUES.F_NAN
      xlon[iday] = !VALUES.F_NAN
      ylat[iday] = !VALUES.F_NAN
   ENDIF ELSE BEGIN
      mld1[iday] = data_insitu[indice_insitu].mean_zhom001
      mld2[iday] = data_insitu[indice_insitu].mean_zhom002
      mld3[iday] = data_insitu[indice_insitu].mean_zhom005
      mld4[iday] = data_insitu[indice_insitu].mean_zhom010
      mld5[iday] = data_insitu[indice_insitu].mean_zhomelch
      std1[iday] = data_insitu[indice_insitu].std_zhom001
      xlon[iday] = data_insitu[indice_insitu].lon
      ylat[iday] = data_insitu[indice_insitu].lat
   ENDELSE
 ENDFOR
;
; compute mean of modeled MLD in boxes around lat,lon insitu data
mymldm = findgen(2,ndays)
mymldr = findgen(2,ndays)
myvarm = findgen(2,ndays)
myvarr = findgen(2,ndays)
;
; define the size of the box
delta=findgen(2)
delta[0] = 0.5
delta[1] = 1.
ndelta=n_elements(delta)
;
; loop on delta
FOR idelta=0, ndelta - 1 DO BEGIN
   ; loop on days
   FOR iday=0, ndays - 1 DO BEGIN
      IF (finite(mld1[iday]) NE 0) THEN BEGIN
         domdef
         mybox = [ -xlon[iday]-delta[idelta] $
                 , -xlon[iday]+delta[idelta] $
                 , ylat[iday]-delta[idelta] $
                 , ylat[iday]+delta[idelta] $
                 ]

         jt = time_insitu[iday]+95
         mldm_box = reform(read_ncdf('mld', jt, /timestep, file=project_id_env+'eddy_mld.nc',  /nostruct, box = mybox))
         mldr_box = reform(read_ncdf('mld', jt, /timestep, file=project_id_env+'noeddy_mld.nc',  /nostruct, box = mybox))

         mymldm[idelta,iday] = moyenne(mldm_box, box = mybox, 'xy')
         mymldr[idelta,iday] = moyenne(mldr_box, box = mybox, 'xy')
         myvarm[idelta,iday] = stddev(mldm_box)
         myvarr[idelta,iday] = stddev(mldr_box)
      ENDIF ELSE BEGIN
         mymldm[idelta,iday] = !VALUES.F_NAN
         mymldr[idelta,iday] = !VALUES.F_NAN
         myvarm[idelta,iday] = !VALUES.F_NAN
         myvarr[idelta,iday] = !VALUES.F_NAN
      ENDELSE
   ENDFOR
 ENDFOR
;
; plot:
; selection de la partie qui nous interesse:
ind = where(finite(mld1) NE 0)
mld1 = mld1[ind]
std1 = std1[ind]
mymldm_1d = mymldm[0, *]
mymldm_1d = mymldm[ind]
mymldr_1d = mymldr[0, *]
mymldr_1d = mymldr[ind]

myvarm_1d = myvarm[0, *]
myvarm_1d = myvarm[ind]
myvarr_1d = myvarr[0, *]
myvarr_1d = myvarr[ind]

time=julday(01,01,2001)+time_insitu
time_save = time
jpt_save = n_elements(time)
time = time[ind]
jpt=n_elements(time)


;++ ok openps, project_od_env+'fig_mld_valid_sta.ps'
plt1d,  -mld1,  't',  -300,  0,  ytitle="meters", title='MLD along cruise track', color = 10, small = [1, 2, 1]
polyfill, [time, reverse(time)], [(-mymldm_1d-myvarm_1d),reverse(-mymldm_1d+myvarm_1d)],  color = 250
polyfill, [time, reverse(time)], [(-mymldr_1d-myvarr_1d),reverse(-mymldr_1d+myvarr_1d)],  color = 150
polyfill, [time, reverse(time)], [(-mld1-std1),reverse(-mld1+std1)], color = 90
plt1d,  -mymldr_1d,  't',  /ov1d, line = 2
plt1d,  -mld1,  't', /ov1d
plt1d,  -mymldm_1d,  't', /ov1d,  line = 1

time = time_save
jpt = jpt_save
plt1d,  ylat,  't',  title='', ytitle='latitude',  color=250, small = [3, 2, 4],  /noerase
plt1d,  xlon,  't',  title='', ytitle='longitude', color=90 , small = [3, 2, 6],  /noerase
closeps

;
; Write results in ${PROJECT_OD}/fig_mld.txt
fullfilename=project_od_env+'fig_mld.txt'
OPENW, lun, fullfilename, /GET_LUN
;++ hanlde error
;
header='jt xlon xlat mld1 mld2 mld3 mld4 mld5 mymldm_0 mymldr_0 myvarm_0 myvarr_0 mymldm_1 mymldr_1 myvarm_1 myvarr_1'
; write header
printf, lun, header
;
; write data lines
FOR iday=0, ndays - 1 DO BEGIN
   printf, lun $
         , iday $
         , xlon[iday] $
         , ylat[iday] $
         , mld1[iday] $
         , mld2[iday] $
         , mld3[iday] $
         , mld4[iday] $
         , mld5[iday] $
         , mymldm[0,iday] $
         , mymldr[0,iday] $
         , myvarm[0,iday] $
         , myvarr[0,iday] $
         , mymldm[1,iday] $
         , mymldr[1,iday] $
         , myvarm[1,iday] $
         , myvarr[0,iday] $
         , mymldm[1,iday] $
         , mymldr[1,iday] $
         , myvarm[1,iday] $
         , myvarr[1,iday]
ENDFOR
;
; close file
FREE_LUN, lun

END