source: trunk/src/fig_mld.pro @ 4

;+
;
;
; .. 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 ++
;
; :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`
;
; :todo:
;
; make it work
;
; consolidate
;
; delta array of 2 : .5 1
;
; noter dans la doc data que l'axe des temps des fichiers modeles est incorrect
;
; :history:
;
; - fplod 20101129T151939Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * 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
      print, '1 ', iday, ' no insitu data'
      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
      print, '2 ', iday, ' 1 insitu data'
      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.
;
; loop on idelta ++
idelta = 0
;
; loop on days
FOR iday=0, ndays - 1 DO BEGIN
      print, ' 0 iday mld1 xlon', iday, mld1[iday], xlon[iday]
   IF (finite(mld1[iday]) NE 0) THEN BEGIN
      print, ' 3 ', iday, ' insitu data'
      print, ' 3a ', 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
      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)
      myvarr[idelta,iday] = stddev(mldr)
   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
;
;write, jt, xlon, ylat, mld1, mld2, mld3, mld4, mld5, mymldm, mymldr,myvarm, myvarr,mymldm[1], mymldr[1],myvarm[1], myvarr[1]
END