Index: /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 496)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 497)
@@ -133,5 +133,5 @@
       real time_step, t_wrt, t_ops
 
-!      REAL rdayvrai,rdaym_ini
+      REAL rdaym_ini
 ! jD_cur: jour julien courant
 ! jH_cur: heure julienne courante
@@ -359,4 +359,6 @@
       CALL geopot  ( ip1jmp1, tsurpk  , pk , pks,  phis  , phi   )
 
+           rdaym_ini  = itau * dtvr / daysec
+
       time = jD_cur + jH_cur
       CALL caldyn 
@@ -435,6 +437,4 @@
          endif
 
-!           rdaym_ini  = itau * dtvr / daysec
-!           rdayvrai   = rdaym_ini  + day_ini
            jD_cur = jD_ref + day_ini - day_ref +                        &
      &          int (itau * dtvr / daysec) 
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 496)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 497)
@@ -146,4 +146,5 @@
       INTEGER :: an, mois, jour
       REAL :: secondes
+      real :: rdaym_ini
 
       LOGICAL first,callinigrads
@@ -626,4 +627,6 @@
 !      CALL FTRACE_REGION_BEGIN("caldyn")
       time = jD_cur + jH_cur 
+           rdaym_ini  = itau * dtvr / daysec
+
 ! ADAPTATION GCM POUR CP(T)
 !      CALL caldyn_p 
Index: /trunk/UTIL/NCL/README
===================================================================
--- /trunk/UTIL/NCL/README	(revision 497)
+++ /trunk/UTIL/NCL/README	(revision 497)
@@ -0,0 +1,32 @@
+This is a set of scripts designed to plot GCM files using NCL graphics.
+=======================================================================
+First version: january 2012, S. Lebonnois
+
+Remarks:
+--------
+
+- Requires NCL to be installed on your system.
+
+- Your simulation file need time, pressure, lon and lat dimensions. 
+The names of the dimensions matter, so you may get an error message if yours are not supported. Just let me know and we'll get this fixed.
+
+- The second variable only gives contours on the 2D maps. 
+Not yet implemented for 1D nor for projections.
+Needs additional options to customize... 
+
+- Hopefully this script will improve with time.
+
+visu-ncl
+--------
+
+- copy the "visu-ncl" bash script to you bin directory 
+
+- modify this copy so that the NCLDIR environment variable correctly point to the UTIL/NCL directory
+
+- when you are in your simulation directory, do :
+   visu-ncl -h |more             : list of available options
+   visu-ncl -f <yoursimu>.nc     : list of variables in <yoursimu>.nc
+
+If you get errors or think of additional features that could be useful, let me know, I'll be happy to improve the script for wider use when possible.
+sllmd@lmd.jussieu.fr
+
Index: /trunk/UTIL/NCL/help.txt
===================================================================
--- /trunk/UTIL/NCL/help.txt	(revision 497)
+++ /trunk/UTIL/NCL/help.txt	(revision 497)
@@ -0,0 +1,112 @@
+************ Usage **************
+
+visu-ncl -f infilename.nc [options] 
+
+ -f infilename.nc   : Input file, mandatory
+
+  This script is designed for longitude/latitude/pressure/time nc files.
+  Required dimension names are: 
+    Time indexes   : Time, time, time_counter
+    Vertical levels: lev, presnivs
+    Latitudes      : lat, latitude
+    Longitudes     : lon, longitudes
+    
+== Options ==
+
+[-h]                : This user manual
+
+[-med media]        : media for output
+
+  Default: x11
+  Possible values:
+     x11        => Display graph on screen
+     ps,eps,pdf => Display graph in a ps,eps or pdf file
+
+[-oname outputname] : name for the output
+
+  Default: quickview (will be the name of the X11 window)
+  The name extension (.ps, .eps, .pdf) is added automatically.
+
+[-type typeplot]    : type of plot
+
+  Default: coupe
+  Possible values:
+     coupe  => Display 1D graph or 2D contour plot
+               You need to give at least 1 varying dimension, 
+               otherwise default give only one point.
+     ortho  => Display 2D orthographic projection
+               You will be defaulted to a full lat/lon map.
+     stereo => Display 2D stereographic projection
+               You will be defaulted to a full lon map.
+
+Options to define the plotted variable:
+---------------------------------------
+
+[-v varname]                     : Variable name for plot
+
+  Default: liste => list the variables in file
+
+[-v2 varname]                    : Variable name to be overplotted
+
+  Default: no second variable
+
+[-tc ntcutmin[,ntcutmax[,tavg]]] : INDEXES for time slice
+
+  Default: ntcutmin=0 / ntcutmax=ntcutmin / tavg=ntcutmax
+  ntcutmin=-999: average over all times
+  tavg=-999    : average between ntcutmin and ntcutmax
+  Known problem: on given files, the default has a problem. 
+  In this case, try a time point different from 0, it should work.
+  
+[-pc pcutmin[,pcutmax[,pavg]]]   : VALUES for pressure slice
+
+  Default: pcutmin=1.e4 / pcutmax=pcutmin / pavg=pcutmax
+  pcutmin=-999 : average over all pressures
+  pavg=-999    : average between pcutmin and pcutmax
+
+[-lc lcutmin[,lcutmax[,lavg]]]   : VALUES for latitudinal slice
+
+  Default: lcutmin=0 / lcutmax=lcutmin / lavg=lcutmax
+  lcutmin=-999 : average over all latitudes
+  lavg=-999    : average between lcutmin and lcutmax
+
+[-Lc Lcutmin[,Lcutmax[,Lavg]]]   : VALUES for longitudinal slice
+
+  Default: Lcutmin=0 / Lcutmax=Lcutmin / Lavg=Lcutmax
+  Lcutmin=-999 : average over all longitudes
+  Lavg=-999    : average between Lcutmin and Lcutmax
+
+Cosmetic options for figure:
+----------------------------
+
+[-lim min,max[,step]] : limits for variable values
+
+  Default: automatic
+  1d     : min and max values for axis (step unused)
+  2d     : min, max and step values for contours
+  If only one value is given, back to default.
+
+Cosmetic options for projections:
+---------------------------------
+
+[-center lon,lat] : longitude and latitude (degrees) of center point
+
+  Default: 
+  For orthographic: 0,40
+  For stereographic: 0,90 (i.e. Northern Hemisphere)
+  To plot Southern Hemisphere, ask for 0,-90 (any other will go back to default)
+  If only one value is given, both are identical.
+  
+[-grid T/F]       : logical to plot the grid or not
+
+  Default: T
+
+[-gcol gridcol]   : color for the grid
+
+  Default: white
+
+[-gspc Dlon,Dlat] : intervals (degrees) for longitude and latitude lines
+
+  Default: 30,20
+  If only one value is given, both are identical.
+
Index: /trunk/UTIL/NCL/planeto.ncl
===================================================================
--- /trunk/UTIL/NCL/planeto.ncl	(revision 497)
+++ /trunk/UTIL/NCL/planeto.ncl	(revision 497)
@@ -0,0 +1,465 @@
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
+load "~/Graphiques/Ncl/visu-utils.ncl"
+
+begin
+
+;---------------------
+;---- load the file
+;---------------------
+
+infile=addfile(inFile,"r")
+
+;---------------------
+;---- get variable list and dimension names
+;---------------------
+
+liste=getfilevarnames(infile)
+nvars=dimsizes(liste)
+
+;---- finding the dimensions 
+tname="dummy"
+pname="dummy"
+lname="dummy"
+Lname="dummy"
+do i=0,nvars-1
+  if ((liste(i).eq."time").or.(liste(i).eq."time_counter").or.(liste(i).eq."Time")) then
+    tname=liste(i)
+  end if
+  if ((liste(i).eq."lev").or.(liste(i).eq."presnivs")) then
+    pname=liste(i)
+  end if
+  if ((liste(i).eq."lat").or.(liste(i).eq."latitude")) then
+    lname=liste(i)
+  end if
+  if ((liste(i).eq."lon").or.(liste(i).eq."longitude")) then
+    Lname=liste(i)
+  end if
+end do
+if ((tname.eq."dummy").or.(pname.eq."dummy").or.(tname.eq."dummy").or.(tname.eq."dummy")) then
+  print("One dimension was not found... See -h for help.")
+  exit
+else
+  names=(/tname,pname,lname,Lname/)
+end if
+
+;---- temporal axis : for Titan, Ls ; for Venus, Vdays.
+Vday=1.0087e7      ; Venus day
+Tday=1.37889e6     ; Titan day
+axels=False
+timeax=infile->$tname$
+timeaxis=tofloat(timeax-timeax(0))/Vday
+copy_VarMeta(timeax,timeaxis)
+timeaxis@units="Local days"
+if (isfilevar(infile,"ls")) then
+   delete(timeaxis)
+   timeaxis=infile->ls
+   axels=True
+end if
+if (isfilevar(infile,"Ls")) then
+   delete(timeaxis)
+   timeaxis=infile->Ls
+   axels=True
+end if
+ndls=dimsizes(dimsizes(timeaxis))
+if (ndls.ne.1) then
+   delete(timeax)
+   timeax=timeaxis(:,0,0)
+   delete(timeaxis)
+   timeaxis=timeax
+end if
+
+;---------------------
+;---- print only list of variables
+;---------------------
+if (var.eq."liste") then
+  print(liste)
+  exit
+end if
+
+;==============================
+; Prepare variable(s) for the plot
+;==============================
+
+;---------------------
+;---- load variable(s)
+;---------------------
+  variable=infile->$var$
+  nbdim=dimsizes(filevardimsizes(infile,var))
+
+  overplot=False
+  diffdims=False
+  if (var2.ne."dummy") then
+     variable2=infile->$var2$
+     nbdim2=dimsizes(filevardimsizes(infile,var2))
+     if (nbdim2.ne.nbdim) then
+       diffdims=True
+     end if
+     overplot=True
+  end if
+  
+;---------------------
+;---- trim variable to desired region
+;---------------------
+     dimname=new(nbdim,string)
+     mindimval=new(nbdim,float)
+     maxdimval=new(nbdim,float)
+     axetim=(/0./)          ; dummy
+     dimavg=new(4,logical)
+     do i=0,3
+       dimavg(i)=False
+     end do
+
+; locate dimensions of the variable
+     flag=positionDims(variable,nbdim,names)
+     if (diffdims) then
+       flag2=positionDims(variable2,nbdim2,names)
+       dimname2=new(nbdim2,string)
+       mindimval2=new(nbdim2,float)
+       maxdimval2=new(nbdim2,float)
+       dimavg2=new(4,logical)
+     end if
+
+; store dim names, associated limits and flags
+     index=0
+     if (flag(0).ne.10) then
+        dimname(index)=tname
+        if (ntcutmin.ne.-999) then
+          mindimval(index)=tofloat(variable&$tname$(ntcutmin))
+          maxdimval(index)=tofloat(variable&$tname$(ntcutmax))
+; A REVOIR !!!
+     ; !! stupid bug, due to precision when time axis in double in file...
+     ;     if ((ntcutmin.eq.0).and.(ntcutmax.eq.0)) then
+     ;       mindimval(index)=mindimval(index)*1.000001
+     ;       maxdimval(index)=maxdimval(index)*1.000001
+     ;     end if
+     ; !! stupid bug, due to precision when time axis in double in file...
+        else
+          mindimval(index)=-999
+          maxdimval(index)=-999
+        end if
+        if (tavg.eq.-999) then
+          dimavg(index)=True
+        else
+     ; trim time axis
+          delete(axetim)
+          axetim=timeaxis({$dimname(0)$|mindimval(0):maxdimval(0)})
+        end if
+        index=index+1
+     end if
+     if (flag(1).ne.10) then
+        dimname(index)=pname
+        mindimval(index)=pcutmin
+        maxdimval(index)=pcutmax
+        if (pavg.eq.-999) then
+          dimavg(index)=True
+        end if
+        index=index+1
+     end if
+     if (flag(2).ne.10) then
+        dimname(index)=lname
+        mindimval(index)=lcutmin
+        maxdimval(index)=lcutmax
+        if (lavg.eq.-999) then
+          dimavg(index)=True
+        end if
+        index=index+1
+     end if
+     if (flag(3).ne.10) then
+        dimname(index)=Lname
+        mindimval(index)=Lcutmin
+        maxdimval(index)=Lcutmax
+        if (Lavg.eq.-999) then
+          dimavg(index)=True
+        end if
+     end if
+
+; if global average, prepare limits
+     do i=0,nbdim-1
+        if (mindimval(i).eq.-999) then
+          mindimval(i)=tofloat(variable&$dimname(i)$(0))
+          maxdimval(i)=tofloat(variable&$dimname(i)$(dimsizes(variable&$dimname(i)$)-1))
+        end if 
+     end do
+
+; if orthographic projection, default limits for lat and lon
+     if (typeplot.eq."ortho") then
+      do i=0,nbdim-1
+        if ((dimname(i).eq.lname).or.(dimname(i).eq.Lname)) then
+          mindimval(i)=tofloat(variable&$dimname(i)$(0))
+          maxdimval(i)=tofloat(variable&$dimname(i)$(dimsizes(variable&$dimname(i)$)-1))
+        end if 
+      end do
+      print("Latitude and Longitude are defaulted to full map.")
+     end if
+     
+; if stereographic projection, default limits for lon
+     if (typeplot.eq."stereo") then
+      do i=0,nbdim-1
+        if (dimname(i).eq.Lname) then
+          mindimval(i)=tofloat(variable&$dimname(i)$(0))
+          maxdimval(i)=tofloat(variable&$dimname(i)$(dimsizes(variable&$dimname(i)$)-1))
+        end if 
+      end do
+      print("Longitude is defaulted to full map.")
+     end if
+     
+; case of second variable
+     if (.not.diffdims) then
+        dimname2  =dimname
+	mindimval2=mindimval
+	maxdimval2=maxdimval
+	dimavg2   =dimavg
+     else
+      index=0
+      if (flag2(0).ne.10) then
+        dimname2(index)=tname
+        index=index+1
+      end if
+      if (flag2(1).ne.10) then
+        dimname2(index)=pname
+        index=index+1
+      end if
+      if (flag2(2).ne.10) then
+        dimname2(index)=lname
+        index=index+1
+      end if
+      if (flag2(3).ne.10) then
+        dimname2(index)=Lname
+      end if
+      do i=0,nbdim2-1
+        do j=0,nbdim-1
+	  if (dimname2(i).eq.dimname(j)) then
+	    mindimval2(i)=mindimval(j)
+	    maxdimval2(i)=maxdimval(j)
+	    dimavg2(i)   =dimavg(j)
+	  end if
+	end do
+      end do
+     end if
+     
+; trim variable(s)
+     vartrim=trimVar(variable,nbdim,dimname,mindimval,maxdimval)
+
+     if (overplot) then
+       vartrim2=trimVar(variable2,nbdim2,dimname2,mindimval2,maxdimval2)
+     end if
+
+;---------------------
+;---- make averages  
+; (if limits are identical, average need to be made to trim the dimension)
+;---------------------
+     varavg=avgVar(vartrim,nbdim,dimavg)
+
+     if (overplot) then
+       varavg2=avgVar(vartrim2,nbdim2,dimavg2)
+     end if
+
+; after trim, redo dim sizes
+     nbdimavg=dimsizes(dimsizes(varavg))
+     if (overplot) then
+       nbdimavg2=dimsizes(dimsizes(varavg2))
+       if (nbdimavg2.ne.nbdimavg) then
+         print("This case is not yet possible: the second variable must be plotted similarly to the first one...")
+         overplot=False
+       end if
+     end if
+
+;==============================
+; Single point
+;==============================
+
+if ((nbdimavg.eq.1).and.(dimsizes(varavg).eq.1)) then
+  print(mindimval)
+  print(varavg)
+  if (overplot) then
+    print(varavg2)
+  end if
+  exit
+end if
+
+;==============================
+;       TRACES
+;==============================
+
+; after trim, redo flag for dims
+     flag=positionDims(varavg,nbdimavg,names)
+
+; 2d maps (e.g. for ortho proj):
+; revert longitudes if decreasing
+  if ( (nbdimavg.eq.2).and.(varavg!0.eq.lname).and.(varavg!1.eq.Lname).and.(varavg&$Lname$(0).gt.varavg&$Lname$(1)) ) then
+     newvaravg=revertLon_2dmap(varavg)
+     delete(varavg)
+     varavg=newvaravg
+     if (overplot) then
+        newvaravg=revertLon_2dmap(varavg2)
+        delete(varavg2)
+        varavg2=newvaravg
+     end if
+  end if
+  
+; Pressure axis used ?
+  axep=False
+  if (flag(1).ne.10) then
+    axep=True
+  end if
+
+; Time axis used ?
+  axet=False
+  if (flag(0).ne.10) then
+    axet=True
+  end if
+
+; limits for contours (2d) or variable axis (1d) ?
+  limits=new(3,float)
+  limits(0)=valmin
+  limits(1)=valmax
+  limits(2)=step
+  if (overplot) then
+    limits2=new(3,float)
+    limits2(0)=valmin2
+    limits2(1)=valmax2
+    limits2(2)=step2
+  end if
+
+; output workstation
+  wks   = gsn_open_wks(media,outputname)
+  resbase          = True 
+  resbase@gsnDraw  = False
+  resbase@gsnFrame = False
+
+;---------------------
+if (typeplot.eq."coupe") then
+;---------------------
+
+; Too many dims
+  if (nbdimavg.gt.2) then
+    print("Too many dimensions for a simple 1D/2D plot...")
+    exit
+  end if
+
+; 2D MAPS
+  if (nbdimavg.eq.2) then
+    res2d = resbase 
+  ;Time axis used with lat or pressure: it must be on X 
+    if ((axet).and.(flag(3).eq.10)) then
+      revvar=inversDim_2d(varavg)
+      delete(varavg)
+      varavg=revvar
+    end if
+  ;defaults and options for plot
+    res2d=deflimits_2d(varavg,res2d)
+    res2d=contours_2d(res2d,limits)
+  ;pressure axis
+    if (axep) then
+      res2d=axePress_2d(res2d)
+    end if
+  ;time axis
+    if ((axet).and.(flag(3).eq.10)) then
+      if (axels) then
+        res2d=axeTimeXls_2d(res2d,axetim)
+      else
+        res2d=axeTimeX_2d(res2d,axetim)
+      end if
+    end if
+    if ((axet).and.(flag(3).ne.10)) then
+      if (axels) then
+        res2d=axeTimeYls_2d(res2d,axetim)
+      else
+        res2d=axeTimeY_2d(res2d,axetim)
+      end if
+    end if
+  ;creating map
+    plotvar=Map_2d(varavg,wks,res2d)
+  ;overlay
+    if (overplot) then
+      res2d_op=res2d
+      res2d_op=contours_2d(res2d_op,limits2)
+      res2d_op=optionsOverPlot_2d(res2d_op)
+      plotvar2=OPlot_2d(varavg2,wks,res2d_op)
+      overlay(plotvar,plotvar2)
+    end if
+  end if
+
+; 1D PLOT
+  if (nbdimavg.eq.1) then
+    res1d = resbase 
+    axename=getvardims(varavg)
+    axe1d=varavg&$axename(0)$
+  ;time axis
+    if (axet) then
+      delete(axe1d)
+      axe1d=axetim
+      if (axels) then
+        res1d=axeTimels_1d(res1d,axe1d)
+      else
+        res1d=axeTime_1d(res1d,axe1d)
+      end if
+    end if
+  ;creating plot
+    plotvar=Plot_1d(varavg,wks,res1d,axe1d,axep,limits)
+  end if
+
+end if
+
+;---------------------
+if (typeplot.eq."ortho") then
+;---------------------
+
+; Need 2D
+  if ((nbdimavg.ne.2).and.((flag(1).eq.10).or.(flag(2).eq.10))) then
+    print("You need 2D lat/lon for orthographic projection...")
+    exit
+  end if
+
+; options
+  center  = (/Lcenter,lcenter/)
+  grd     = grid
+  spacing = (/gLspc,glspc/)
+  gridcol = gcol
+  
+  res2d = resbase 
+  res2d = contours_2d(res2d,limits)
+  res2d = optionsOrtho(res2d,center,grd,spacing,gridcol)
+
+;creating plot
+  plotvar=Ortho_2d(varavg,wks,res2d)
+
+end if
+
+;---------------------
+if (typeplot.eq."stereo") then
+;---------------------
+
+; Need 2D
+  if ((nbdimavg.ne.2).and.((flag(1).eq.10).or.(flag(2).eq.10))) then
+    print("You need 2D lat/lon for stereographic projection...")
+    exit
+  end if
+
+; options
+  if (lcenter.eq.-90) then
+    hemisph="SH"
+  else
+    hemisph="NH"
+  end if
+  grd     = grid
+  spacing = (/gLspc,glspc/)
+  gridcol = gcol
+  
+  res2d = resbase 
+  res2d = contours_2d(res2d,limits)
+  res2d = optionsStereo(res2d,hemisph,grd,spacing,gridcol)
+  
+;creating plot
+  plotvar=Stereo_2d(varavg,wks,res2d)
+
+end if
+
+;---------------------
+  draw(plotvar)                   
+  frame(wks)
+
+end
+
Index: /trunk/UTIL/NCL/visu-ncl
===================================================================
--- /trunk/UTIL/NCL/visu-ncl	(revision 497)
+++ /trunk/UTIL/NCL/visu-ncl	(revision 497)
@@ -0,0 +1,221 @@
+#!/bin/bash
+
+# Directory where the NCL files are located ** EDIT **
+
+export NCLDIR=~/LMDZ5/UTIL/NCL
+
+########################################################################
+# Mise en place des valeurs par defaut
+########################################################################
+
+inFile="nofile"
+var="liste"
+var2="dummy"
+media="x11"
+outputname="quickview"
+typeplot="coupe"
+ntcutmin=0
+ntcutmax=$ntcutmin
+tavg=-999
+pcutmin=1.e4
+pcutmax=$pcutmin
+pavg=-999
+lcutmin=0
+lcutmax=$lcutmin
+lavg=-999
+Lcutmin=0
+Lcutmax=$Lcutmin
+Lavg=-999
+# (if limits are identical, average need to be made to trim the dimension)
+
+# -888 => default variable chosen by NCL
+valmin=-888
+valmax=-888
+step=-888
+valmin2=-888
+valmax2=-888
+step2=-888
+
+Lcenter=0
+lcenter=40
+grid="True"
+gLspc=30
+glspc=20
+gcol="white"
+
+########################################################################
+# lecture des options
+########################################################################
+
+while (($# > 0))
+  do
+  case $1 in
+      "-h") cat $NCLDIR/help.txt
+	  exit;;
+      "-f")
+	  inFile=$2 ; shift ; shift ;;
+      "-v")
+	  var=$2 ; shift ; shift ;;
+      "-v2")
+	  var2=$2 ; shift ; shift ;;
+      "-med")
+	  media=$2 ; shift ; shift ;;
+      "-oname")
+	  outputname=$2 ; shift ; shift ;;
+      "-type")
+          typeplot=$2
+          if [ $typeplot = "stereo" ]; then
+            lcenter=90
+          fi
+	  shift ; shift ;;
+      "-tc")
+	  ntcutmin=${2%%,*}
+          tmp=${2#*,}
+          ntcutmax=${tmp%%,*}
+          if [ $ntcutmin = $ntcutmax ]; then
+            tavg=-999
+          else
+            tavg=${tmp#*,}
+          fi
+          shift ; shift ;;
+      "-pc")
+	  pcutmin=${2%%,*}
+          tmp=${2#*,}
+          pcutmax=${tmp%%,*}
+          if [ $pcutmin = $pcutmax ]; then
+            pavg=-999
+          else
+            pavg=${tmp#*,}
+          fi
+          shift ; shift ;;
+      "-lc")
+	  lcutmin=${2%%,*}
+          tmp=${2#*,}
+          lcutmax=${tmp%%,*}
+          if [ $lcutmin = $lcutmax ]; then
+            lavg=-999
+          else
+            lavg=${tmp#*,}
+          fi
+          shift ; shift ;;
+      "-Lc")
+	  Lcutmin=${2%%,*}
+          tmp=${2#*,}
+          Lcutmax=${tmp%%,*}
+          if [ $Lcutmin = $Lcutmax ]; then
+            Lavg=-999
+          else
+            Lavg=${tmp#*,}
+          fi
+          shift ; shift ;;
+       "-lim")
+	  valmin=${2%%,*}
+          tmp=${2#*,}
+          valmax=${tmp%%,*}
+          if [ $valmin = $valmax ]; then
+          # si pas les deux, remises au default
+             valmin=-888
+             valmax=-888
+          else
+            step=${tmp#*,}
+          fi
+          if [ $step = $valmax ]; then
+          # si pas de step, il est remis au default
+             step=-888
+          fi
+          shift ; shift ;;
+       "-lim2")
+	  valmin2=${2%%,*}
+          tmp=${2#*,}
+          valmax2=${tmp%%,*}
+          if [ $valmin2 = $valmax2 ]; then
+          # si pas les deux, remises au default
+             valmin2=-888
+             valmax2=-888
+          else
+            step2=${tmp#*,}
+          fi
+          if [ $step2 = $valmax2 ]; then
+          # si pas de step, il est remis au default
+             step2=-888
+          fi
+          shift ; shift ;;
+       "-center")
+	  Lcenter=${2%%,*}
+          tmp=${2#*,}
+          lcenter=${tmp%%,*}
+          shift ; shift ;;
+       "-grid")
+	  if [ $2 = "T" -o $2 = "y" -o $2 = "Y" ]; then
+	    grid="True"
+	  fi
+	  if [ $2 = "F" -o $2 = "n" -o $2 = "N" ]; then
+	    grid="False"
+	  fi
+	  shift ; shift ;;
+      "-gcol")
+	  gcol=$2 ; shift ; shift ;;
+       "-gspc")
+	  gLspc=${2%%,*}
+          tmp=${2#*,}
+          glspc=${tmp%%,*}
+	  shift ; shift ;;
+  esac
+done
+
+###############################################################
+if [ $inFile = "nofile" ]; then
+	echo "You have to specify a file ! ( -f filename )"; exit
+fi 
+
+if [ -f $inFile ]; then
+	echo "File "$inFile" exists"
+else 
+	echo "File "$inFile" does not exists"; exit
+fi 
+
+if [ $typeplot = "ortho" ]; then
+    # lat and lon put to 2d, defaulting to full map in planeto.ncl
+	lcutmin=-90
+        lcutmax=90
+	lavg=0
+        Lcutmin=0
+        Lcutmax=180
+	Lavg=0
+fi 
+
+export SIMUDIR=`pwd`
+inFile=$SIMUDIR/$inFile
+cd $NCLDIR
+ncl inFile=\"$inFile\"  \
+    var=\"$var\" \
+    var2=\"$var2\" \
+    media=\"$media\" \
+    outputname=\"$outputname\" \
+    typeplot=\"$typeplot\" \
+    ntcutmin=$ntcutmin \
+    ntcutmax=$ntcutmax \
+    tavg=$tavg \
+    pcutmin=$pcutmin \
+    pcutmax=$pcutmax \
+    pavg=$pavg \
+    lcutmin=$lcutmin \
+    lcutmax=$lcutmax \
+    lavg=$lavg \
+    Lcutmin=$Lcutmin \
+    Lcutmax=$Lcutmax \
+    Lavg=$Lavg \
+    valmin=$valmin \
+    valmax=$valmax \
+    step=$step \
+    valmin2=$valmin2 \
+    valmax2=$valmax2 \
+    step2=$step2 \
+    Lcenter=$Lcenter \
+    lcenter=$lcenter \
+    grid=$grid \
+    gLspc=$gLspc \
+    glspc=$glspc \
+    gcol=\"$gcol\" \
+    planeto.ncl
+
Index: /trunk/UTIL/NCL/visu-utils.ncl
===================================================================
--- /trunk/UTIL/NCL/visu-utils.ncl	(revision 497)
+++ /trunk/UTIL/NCL/visu-utils.ncl	(revision 497)
@@ -0,0 +1,629 @@
+;;;;;;;;;;;; positionDims ;;;;;;;;;;;;;;;;
+
+undef("positionDims")
+function positionDims(var,nbdim,names)
+;============================
+
+; var   : variable
+; nbdim : number of dimensions for var
+; names : names of the 4 dimensions: time,levs,lat,lon
+
+begin
+
+  flag=(/10,10,10,10/)
+  do i=0,nbdim-1
+     if (var!i.eq.names(0)) then
+        flag(0)=i
+     end if
+     if (var!i.eq.names(1)) then
+        flag(1)=i
+     end if
+     if (var!i.eq.names(2)) then
+        flag(2)=i
+     end if
+     if (var!i.eq.names(3)) then
+        flag(3)=i
+     end if
+  end do
+
+  return(flag)
+
+end
+
+;;;;;;;;;;;; trimVar ;;;;;;;;;;;;;;;;
+
+undef("trimVar")
+function trimVar(var,nbdim,dimname,mindimval,maxdimval)
+;============================
+
+; var      : variable
+; nbdim    : number of dimensions for var
+; dimname  : names of the dimensions for this var
+; mindimval: min values for the dimensions
+; maxdimval: max values for the dimensions
+
+begin
+
+  if (nbdim.eq.1) then
+     vartrim=var({$dimname(0)$|mindimval(0):maxdimval(0)})
+  end if
+  if (nbdim.eq.2) then
+       vartrim=var({$dimname(0)$|mindimval(0):maxdimval(0)},{$dimname(1)$|mindimval(1):maxdimval(1)})
+  end if
+  if (nbdim.eq.3) then
+       vartrim=var({$dimname(0)$|mindimval(0):maxdimval(0)},{$dimname(1)$|mindimval(1):maxdimval(1)},{$dimname(2)$|mindimval(2):maxdimval(2)})
+  end if
+  if (nbdim.eq.4) then
+       vartrim=var({$dimname(0)$|mindimval(0):maxdimval(0)},{$dimname(1)$|mindimval(1):maxdimval(1)},{$dimname(2)$|mindimval(2):maxdimval(2)},{$dimname(3)$|mindimval(3):maxdimval(3)})
+  end if
+
+  return(vartrim)
+
+end
+
+;;;;;;;;;;;; avgVar ;;;;;;;;;;;;;;;;
+
+undef("avgVar")
+function avgVar(vartrim,nbdim,dimavg)
+;============================
+
+; vartrim : variable
+; nbdim   : number of dimensions for var
+; dimavg  : logicals for averaging the dimensions
+
+begin
+
+  if ((nbdim.eq.4).and.(dimavg(3))) then
+       vartrim_3=dim_avg_n_Wrap(vartrim,3)
+  else
+       vartrim_3=vartrim
+  end if
+  if ((nbdim.ge.3).and.(dimavg(2))) then
+       vartrim_2=dim_avg_n_Wrap(vartrim_3,2)
+  else
+       vartrim_2=vartrim_3
+  end if
+  if ((nbdim.ge.2).and.(dimavg(1))) then
+       vartrim_1=dim_avg_n_Wrap(vartrim_2,1)
+  else
+       vartrim_1=vartrim_2
+  end if
+  if ((nbdim.ge.1).and.(dimavg(0))) then
+       varavg=dim_avg_n_Wrap(vartrim_1,0)
+  else
+       varavg=vartrim_1
+  end if
+
+  return(varavg)
+
+end
+
+;;;;;;;;;;;; deflimits_2d ;;;;;;;;;;;;;;;;
+
+undef("deflimits_2d")
+function deflimits_2d(var,res)
+;============================
+
+; var : 2D variable
+; res : ressource
+
+begin
+
+  dims=dimsizes(var)
+  dimname=getvardims(var)
+  res@trXMinF = tofloat(var&$dimname(1)$(0))
+  res@trXMaxF = tofloat(var&$dimname(1)$(dims(1)-1))
+  res@trYMinF = tofloat(var&$dimname(0)$(0))
+  res@trYMaxF = tofloat(var&$dimname(0)$(dims(0)-1))
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; contours_2d ;;;;;;;;;;;;;;;;
+
+undef("contours_2d")
+function contours_2d(res,limits)
+;============================
+
+; res    : modified ressource for the plot
+; limits : min, max and values for contours (-888 => automatic)
+
+begin
+
+  if (limits(0).ne.-888) then
+     res@cnLevelSelectionMode = "ManualLevels"
+     res@cnMinLevelValF       = limits(0)
+     res@cnMaxLevelValF       = limits(1)
+     if (limits(2).ne.-888) then
+        res@cnLevelSpacingF   = limits(2)
+     end if
+  end if
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; optionsOrtho ;;;;;;;;;;;;;;;;
+
+undef("optionsOrtho")
+function optionsOrtho(res,center,grid,spacing,gridcol)
+;============================
+
+; res    : modified ressource for the plot
+; center : lon,lat coordinates for center of projection
+; grid   : logical for printing grid
+; spacing: lon,lat spacing for grid
+; gridcol: color for grid
+
+begin
+
+  res@mpCenterLonF      = center(0)  ; choose center lon
+  res@mpCenterLatF      = center(1)  ; choose center lat
+  res@mpGridAndLimbOn   = grid
+  res@mpGridLonSpacingF = spacing(0)
+  res@mpGridLatSpacingF = spacing(1)
+  res@mpGridLineColor   = gridcol
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; optionsStereo ;;;;;;;;;;;;;;;;
+
+undef("optionsStereo")
+function optionsStereo(res,hm,grid,spacing,gridcol)
+;============================
+
+; res    : modified ressource for the plot
+; hm     : hemisphere to be plotted
+; grid   : logical for printing grid
+; spacing: lon,lat spacing for grid
+; gridcol: color for grid
+
+begin
+
+  res@gsnPolar          = hm
+
+  res@mpGridAndLimbOn   = grid
+  res@mpGridLonSpacingF = spacing(0)
+  res@mpGridLatSpacingF = spacing(1)
+  res@mpGridLineColor   = gridcol
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; optionsOverPlot_2d ;;;;;;;;;;;;;;;;
+
+undef("optionsOverPlot_2d")
+function optionsOverPlot_2d(res)
+;============================
+
+; res    : modified ressource for the plot
+
+begin
+
+  res@cnLineLabelFormat            = "@5.1f"
+  res@cnLineLabelsOn               = True
+  res@cnLabelMasking               = True		
+  res@cnLineLabelBackgroundColor   = -1
+  res@cnInfoLabelOn                = False
+  res@gsnContourNegLineDashPattern = 1 
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; axePress_2d ;;;;;;;;;;;;;;;;
+
+undef("axePress_2d")
+function axePress_2d(res)
+;============================
+
+; res : modified ressource for the plot
+
+begin
+
+  res@tiYAxisString  = "Pressure"
+  res@trYReverse     = True
+  res@gsnYAxisIrregular2Log = True
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; axePress_1d ;;;;;;;;;;;;;;;;
+
+undef("axePress_1d")
+function axePress_1d(res,axe,limits)
+;============================
+
+; res    : modified ressource for the plot
+; axe    : Y axis for plot (pressure)
+; limits : min and max values for variable axis
+
+begin
+
+  res@tiYAxisString  = "Pressure"
+  res@trYReverse     = True
+  res@xyYStyle       = "Log"
+  res@trYMinF         = axe(0)
+  res@trYMaxF         = axe(dimsizes(axe)-1)
+  if (limits(0).ne.-888) then
+    res@trXMinF       = limits(0)
+    res@trXMaxF       = limits(1)
+  end if
+
+  return(res)
+
+end
+
+
+undef("axeNoPress_1d")
+function axeNoPress_1d(res,axe,limits)
+;============================
+
+; res    : modified ressource for the plot
+; axe    : X axis for plot (not pressure)
+; limits : min and max values for variable axis
+
+begin
+
+  res@trXMinF         = axe(0)
+  res@trXMaxF         = axe(dimsizes(axe)-1)
+  if (limits(0).ne.-888) then
+    res@trYMinF       = limits(0)
+    res@trYMaxF       = limits(1)
+  end if
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; axeTimeX_2d ;;;;;;;;;;;;;;;;
+
+undef("axeTimeX_2d")
+function axeTimeX_2d(res,axe)
+;============================
+
+; res : modified ressource for the plot
+; axe : X axis is time
+
+begin
+
+  res@tiXAxisString        = "Time (local days)"
+  res@sfXArray             = axe
+  res@tmXBLabelFontHeightF = 0.015
+  res@trXMinF              = axe(0)
+  res@trXMaxF              = axe(dimsizes(axe)-1)
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; axeTimeXls_2d ;;;;;;;;;;;;;;;;
+
+undef("axeTimeXls_2d")
+function axeTimeXls_2d(res,axe)
+;============================
+
+; res : modified ressource for the plot
+; axe : X axis is Ls
+
+begin
+
+  res@tiXAxisString            = "Solar longitude"
+  res@sfXArray                 = axe
+  res@gsnXAxisIrregular2Linear = True
+  res@tmXBMode                 = "Manual"	
+  res@tmXBTickStartF           = 0
+  res@tmXBTickEndF             = 360
+  res@tmXBTickSpacingF         = 30
+  res@tmXBMinorPerMajor        = 2
+  res@tmXBLabelFontHeightF     = 0.015
+  res@trXMinF                  = axe(0)
+  res@trXMaxF                  = axe(dimsizes(axe)-1)
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; axeTimeY_2d ;;;;;;;;;;;;;;;;
+
+undef("axeTimeY_2d")
+function axeTimeY_2d(res,axe)
+;============================
+
+; res : modified ressource for the plot
+; axe : Y axis is time
+
+begin
+
+  res@tiYAxisString        = "Time (local days)"
+  res@sfYArray             = axe
+  res@tmYLLabelFontHeightF = 0.015
+  res@trYMinF              = axe(0)
+  res@trYMaxF              = axe(dimsizes(axe)-1)
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; axeTimeYls_2d ;;;;;;;;;;;;;;;;
+
+undef("axeTimeYls_2d")
+function axeTimeYls_2d(res,axe)
+;============================
+
+; res : modified ressource for the plot
+; axe : Y axis is Ls
+
+begin
+
+  res@tiYAxisString            = "Solar longitude"
+  res@sfYArray                 = axe
+  res@gsnYAxisIrregular2Linear = True
+  res@tmYLMode                 = "Manual"	
+  res@tmYLTickStartF           = 0
+  res@tmYLTickEndF             = 360
+  res@tmYLTickSpacingF         = 30
+  res@tmYLMinorPerMajor        = 2
+  res@tmYLLabelFontHeightF     = 0.015
+  res@trYMinF                  = axe(0)
+  res@trYMaxF                  = axe(dimsizes(axe)-1)
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; axeTime_1d ;;;;;;;;;;;;;;;;
+
+undef("axeTime_1d")
+function axeTime_1d(res,axe)
+;============================
+
+; res : modified ressource for the plot
+; axe : X axis is time
+
+begin
+
+  res@tiXAxisString        = "Time (local days)"
+  res@tmXBLabelFontHeightF = 0.015
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; axeTimels_1d ;;;;;;;;;;;;;;;;
+
+undef("axeTimels_1d")
+function axeTimels_1d(res,axe)
+;============================
+
+; res : modified ressource for the plot
+; axe : X axis is Ls
+
+begin
+
+  res@tiXAxisString            = "Solar longitude"
+  res@tmXBMode                 = "Manual"	
+  res@tmXBTickStartF           = 0
+  res@tmXBTickEndF             = 360
+  res@tmXBTickSpacingF         = 30
+  res@tmXBMinorPerMajor        = 2
+  res@tmXBLabelFontHeightF     = 0.015
+
+  return(res)
+
+end
+
+;;;;;;;;;;;; inversDim_2d ;;;;;;;;;;;;;;;;
+
+undef("inversDim_2d")
+function inversDim_2d(var)
+;============================
+
+; var is a 2d variable
+
+begin
+
+  dimname=getvardims(var)
+  tmpvar = var($dimname(1)$|:,$dimname(0)$|:)
+
+  return(tmpvar)
+
+end
+
+;;;;;;;;;;;; revertLon_2dmap ;;;;;;;;;;;;;;;;
+
+undef("revertLon_2dmap")
+function revertLon_2dmap(var)
+;============================
+
+; var is a 2d lat-lon variable
+
+begin
+
+  tmpvar  = var
+  dimname = getvardims(var)
+  Lname   = dimname(1)
+  nbL     = dimsizes(var&$Lname$)
+  tmpaxL  = new(nbL,float)
+
+;revert lon axis:
+  do i=0,nbL-1
+    tmpaxL(i)=var&$Lname$(nbL-1-i)
+  end do
+  tmpvar&$Lname$ = tmpaxL
+
+;revert values along lon axis:
+  do i=0,nbL-1
+    tmpvar(:,i) = var(:,nbL-1-i)
+  end do
+
+  return(tmpvar)
+
+end
+
+;;;;;;;;;;;; Map_2d ;;;;;;;;;;;;;;;;
+
+undef("Map_2d")
+function Map_2d(var,wks,res)
+;============================
+
+; var : 2D variable
+; wks : workstation
+; res : ressource
+
+begin
+
+; main plot characteristics
+
+  res@gsnSpreadColors = True
+  res@tmXTOn          = False
+  res@cnFillOn        = True
+  res@cnLineLabelsOn  = True
+  res@cnLabelMasking  = True
+  res@cnLineLabelBackgroundColor = "white"
+  res@lbOrientation   = "vertical"
+
+; Create plot
+
+  plotvar  = gsn_csm_contour(wks, var, res )
+
+  return(plotvar)
+
+end
+
+;;;;;;;;;;;; OPlot_2d ;;;;;;;;;;;;;;;;
+
+undef("OPlot_2d")
+function OPlot_2d(var,wks,res)
+;============================
+
+; var : 2D variable
+; wks : workstation
+; res : ressource
+
+begin
+
+; main plot characteristics
+
+  res@tmXTOn          = False
+  res@tmYROn          = False
+  res@cnFillOn        = False
+
+  res@gsnLeftString   = ""
+  res@gsnRightString  = ""
+
+; Create plot
+
+  plotvar  = gsn_csm_contour(wks, var, res )
+
+  return(plotvar)
+
+end
+
+;;;;;;;;;;;; Plot_1d ;;;;;;;;;;;;;;;;
+
+undef("Plot_1d")
+function Plot_1d(var,wks,res,axe,axep,limits)
+;============================
+
+; var            : 1D variable
+; wks            : workstation
+; res            : ressource
+; axe            : axis for plot
+; axep           : logical indication for pressure axis
+; limits         : min and max values for variable axis
+
+begin
+
+; main plot characteristics
+
+; Create plot
+; pressure axis ?
+  if (axep) then
+    res=axePress_1d(res,axe,limits)
+    plotvar  = gsn_csm_xy(wks, var, axe, res )
+  else
+    res1d=axeNoPress_1d(res,axe,limits)
+    plotvar  = gsn_csm_xy(wks, axe, var, res )
+  end if
+
+  return(plotvar)
+
+end
+
+;;;;;;;;;;;; Ortho_2d ;;;;;;;;;;;;;;;;
+
+undef("Ortho_2d")
+function Ortho_2d(var,wks,res)
+;============================
+
+; var : 2D variable
+; wks : workstation
+; res : ressource
+
+begin
+
+; main plot characteristics
+
+  res@gsnSpreadColors = True
+  res@tmXTOn          = False
+  res@cnFillOn        = True
+  res@cnLineLabelsOn  = True
+  res@cnLabelMasking  = True
+  res@cnLineLabelBackgroundColor = "white"
+  res@lbOrientation   = "vertical"
+
+  res@mpProjection      = "Orthographic"       ; choose projection
+  res@mpOutlineOn       = False
+  res@mpPerimOn         = False             ; turn off box around plot
+  res@mpFillOn          = False
+  res@mpGridLineDashPattern = 2
+  
+; Create plot
+
+  plotvar  = gsn_csm_contour_map(wks, var, res )
+
+  return(plotvar)
+
+end
+
+;;;;;;;;;;;; Stereo_2d ;;;;;;;;;;;;;;;;
+
+undef("Stereo_2d")
+function Stereo_2d(var,wks,res)
+;============================
+
+; var : 2D variable
+; wks : workstation
+; res : ressource
+
+begin
+
+; main plot characteristics
+
+  res@gsnSpreadColors = True
+  res@tmXTOn          = False
+  res@cnFillOn        = True
+  res@cnLineLabelsOn  = True
+  res@cnLabelMasking  = True
+  res@cnLineLabelBackgroundColor = "white"
+  res@lbOrientation   = "vertical"
+
+  res@mpProjection      = "Stereographic"       ; choose projection
+  res@mpOutlineOn       = False
+  res@mpPerimOn         = False             ; turn off box around plot
+  res@mpFillOn          = False
+  
+; Create plot
+
+  plotvar  = gsn_csm_contour_map_polar(wks, var, res )
+
+  return(plotvar)
+
+end
+
