source: trunk/src/fig_mld.pro @ 5

;+
;
;
; .. program:: fig_mld
;
; ===========
; fig_mld.pro
; ===========
;
; main POMME project program
;
; - read the 2 mld models :file:`${POMME_ID}/noeddy_mld.nc` and :file:`${POMME_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/4° around insitu measurement
; location
; - plot ++
; - write file:`${POMME_OD}/fig_mld.txt`
;
; :examples:
;
; ::
;
;   IDL> fig_mld
;   IDL> exit
;
; :uses:
;
; :func:`read_ncdf`
; :func:`julday`
; :func:`pltt`
; :func:`moyenne`
; :func:`domdef`
;
; :ref:`init_pom.pro`
; :func:`file_prn_to_mem`
; :func:`file_asc_to_mem`
; :func:`dmean_mld`
; :func:`join_dmean_asc`
;
; :func:`reform`
;
; :todo:
;
; make it work
;
; consolidate
;
; better organisation of output files
;
; noter dans la doc data que l'axe des temps des fichiers modeles est incorrect
;
; :history:
;
; - 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
;
@init_pom
;
; read MLD in model files
mldm = -read_ncdf('mld', 0, 383, /timestep, file=iodir+'eddy_mld.nc',  /nostruct)
mldr = -read_ncdf('mld', 0, 383, /timestep, file=iodir+'noeddy_mld.nc',  /nostruct)
;
time=julday(09,27,2000)+dindgen(jpt)
;
; plot modeled MLD
;openps, 'fig_mld.ps' ++
pltt, mldm , box=[-22,-16, 38, 44],  't',  -200,  0, color = 70,  title = 'MLD'
pltt, mldr , box=[-22,-16, 38, 44],  't',  color = 220,  /ov1d
;closeps
;
; 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)
;
; calcul du temps des données insitu avec comme 0 le 20000927 ++
time_insitu=julday(1,1,2001) - julday(09,27,2000) - 1 + data_insitu.doy
;
; build compatible insitu arrays with models arrays
ndays = n_elements(time)
mld1 =  findgen(ndays)
mld2 = findgen(ndays)
mld3 = findgen(ndays)
mld4 = findgen(ndays)
mld5 = findgen(ndays)
xlon = findgen(ndays)
ylat = findgen(ndays)
FOR iday=0, ndays - 1 DO BEGIN
   indice_insitu=WHERE(time_insitu 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
      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
      xlon[iday] = data_insitu[indice_insitu].lon
      ylat[iday] = data_insitu[indice_insitu].lat
   ENDELSE
ENDFOR
;
; plot insitu daily mean MLD ++
;
; compute mean of modeled MLD in box 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
         print, ' 3 iday xlon', iday, xlon[iday]
         domdef
         mybox = [ xlon[iday]-delta[idelta] $
                 , xlon[iday]+delta[idelta] $
                 , ylat[iday]-delta[idelta] $
                 , ylat[iday]+delta[idelta] $
                 ]
         help, mybox
         print, mybox

         ; ++ verifier que la moyenne est bien faite dans la boite et pas dans tout le domaine
         ;++ pb t point
         ;mldm_box = reform(read_ncdf('mld', iday, /timestep, file=iodir+'eddy_mld.nc',  /nostruct, box = mybox))
         ;mldr_box = reform(read_ncdf('mld', iday, /timestep, file=iodir+'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
         ; print, '4 ',iday, ' no insitu data'
         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
;
; Write results in ${POMME_OD}/fig_mld.txt
filename=pomme_od_env+'fig_mld.txt'
OPENW,lun, filename, /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[1,iday]
ENDFOR
;
; close file
FREE_LUN, lun
;
END