Index: trunk/000-USERS
===================================================================
--- trunk/000-USERS	(revision 85)
+++ trunk/000-USERS	(revision 86)
@@ -1,4 +1,2 @@
-Last update: 20/02/2011
-
 ------------------------------------------------------------------------------------------------------------------------
 CONTRIBUTORS [are allowed to commit]
Index: trunk/mars/libf/phymars/dustlift.F
===================================================================
--- trunk/mars/libf/phymars/dustlift.F	(revision 85)
+++ trunk/mars/libf/phymars/dustlift.F	(revision 86)
@@ -1,3 +1,4 @@
-      SUBROUTINE dustlift(ngrid,nlay,nq,rho,pcdh_true,pcdh,co2ice,
+      SUBROUTINE dustlift(ngrid,nlay,nq,rho,
+     $                  pcdh_true,pcdh,co2ice,
      $                  dqslift)
       IMPLICIT NONE
@@ -44,4 +45,21 @@
       DATA stress_seuil/0.0225/   ! stress seuil soulevement (N.m2)
 
+#ifdef MESOSCALE
+!!!! AS: In the mesoscale model we'd like to easily set
+!!!! AS: ... stress for lifting
+!!!! AS: you have to compile with -DMESOSCALE to do so
+      REAL alpha
+      INTEGER ierr
+        OPEN(99,file='stress.def',status='old',form='formatted'
+     .   ,iostat=ierr)
+        !!! no file => default values
+        IF(ierr.EQ.0) THEN
+          READ(99,*) stress_seuil
+          READ(99,*) alpha
+          write(*,*) 'USER-DEFINED threshold: ', stress_seuil, alpha
+          CLOSE(99)
+          alpha_lift(1:nq) = alpha
+        ENDIF
+#endif
 
 c     ---------------------------------
Index: trunk/mars/libf/phymars/meso_dustlift.F
===================================================================
--- trunk/mars/libf/phymars/meso_dustlift.F	(revision 85)
+++ 	(revision )
@@ -1,134 +1,0 @@
-      SUBROUTINE dustlift(ngrid,nlay,nq,rho,
-     $                  pcdh_true,pcdh,co2ice,
-     $                  dqslift)
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c  Dust lifting by surface winds
-c  Computing flux to the middle of the first layer
-c  (Called by vdifc)
-c
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "tracer.h"
-
-c
-c   arguments:
-c   ----------
-
-c   INPUT
-      integer ngrid, nlay, nq  
-      real rho(ngrid)  ! density (kg.m-3) at surface
-      real pcdh_true(ngrid) ! Cd 
-      real pcdh(ngrid) ! Cd * |V| 
-      real co2ice(ngrid)
-
-c   OUTPUT
-      real dqslift(ngrid,nq) !surface dust flux to mid-layer (<0 when lifing)
-c     real pb(ngrid,nlay) ! diffusion to surface coeff.
-
-c   local:
-c   ------
-      INTEGER ig,iq
-      REAL fhoriz(ngridmx)  ! Horizontal dust flux 
-      REAL ust,us
-      REAL stress_seuil
-      SAVE stress_seuil
-c      DATA stress_seuil/0.0225/   ! stress seuil soulevement (N.m2)
-!****WRF
-!****WRF: additional ASCII file to define dust opacity
-          REAL alpha
-          INTEGER ierr
-          OPEN(99,file='stress.def',status='old',form='formatted'
-     .     ,iostat=ierr)
-          IF(ierr.NE.0) THEN
-             stress_seuil = 0.0225
-!!!! defini dans initracer
-!             alpha = 1.
-             write(*,*) 'No file stress.def - set ', 
-     .                     stress_seuil, alpha_lift
-             !stop
-          ELSE
-             READ(99,*) stress_seuil
-             READ(99,*) alpha
-             write(*,*) 'definir seuil stress : ', stress_seuil, alpha
-             CLOSE(99)
-             alpha_lift(1) = alpha
-          ENDIF
-!****WRF
-!****WRF
-
-
-
-c     ---------------------------------
-c     Computing horizontal flux: fhoriz
-c     ---------------------------------
-
-      do ig=1,ngrid
-          fhoriz(ig) = 0.      ! initialisation
-
-c         Selection of points where surface dust is available
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c         if (latid(ig).ge.80.) goto 99  ! N permanent  polar caps
-c         if (latid(ig).le.-80.) goto 99 ! S polar deposits
-c         if  ((longd(ig).ge.-141. .and. longd(ig).le.-127.)
-c    &   .and.(latid(ig).ge.12.   .and. latid(ig).le.23.))goto 99 ! olympus
-c         if  ((longd(ig).ge.-125. .and. longd(ig).le.-118.)
-c    &   .and.(latid(ig).ge.-12.   .and. latid(ig).le.-6.))goto 99 ! Arsia
-c         if  ((longd(ig).ge.-116. .and. longd(ig).le.-109.)
-c    &   .and.(latid(ig).ge.-5.   .and. latid(ig).le. 5.))goto 99 ! pavonis
-c         if  ((longd(ig).ge.-109. .and. longd(ig).le.-100.)
-c    &   .and.(latid(ig).ge. 7.   .and. latid(ig).le. 16.))goto 99 ! ascraeus
-c         if  ((longd(ig).ge.  61. .and. longd(ig).le.  63.)
-c    &   .and.(latid(ig).ge. 63. .and. latid(ig).le. 64.))goto 99 !weird point
-          if (co2ice(ig).gt.0.) goto 99
-
-
-c         Is the wind strong enough ?
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          ust = sqrt(stress_seuil/rho(ig))
-          us = pcdh(ig) /  sqrt(pcdh_true(ig)) ! ustar=cd*v /sqrt(cd)
-          if (us.gt.ust) then
-c            If lifting ? 
-c            Calcul du flux suivant Marticorena ( en fait white (1979))
-
-             fhoriz(ig) = 2.61*(rho(ig)/g) *
-     &      (us -ust) * (us + ust)**2
-          end if
- 99      continue
-      end do
-
-c     -------------------------------------
-c     Computing vertical flux and diffusion
-c     -------------------------------------
- 
-       do iq=1,nq
-         do ig=1,ngrid
-             dqslift(ig,iq)= -alpha_lift(iq)* fhoriz(ig)
-
-
-cc  le  flux vertical remplace le terme de diffusion turb. qui est mis a zero
-c            zb(ig,1) = 0.
-cc           If surface deposition by turbulence diffusion (impaction...)
-cc           if(fhoriz(ig).ne.0) then
-cc           zb(ig,1) = zcdh(ig)*zb0(ig,1)
-cc           AMount of Surface deposition !
-cc           pdqs_dif(ig,iq)=pdqs_dif(ig,iq) +
-cc    &      zb(ig,1)*zq(ig,1,iq)/ptimestep
-cc          write(*,*) 'zb(1)  = ' ,  zb(ig,1),zcdh(ig),zb0(ig,1)
-cc
-
-         enddo
-       enddo
-
-      RETURN
-      END
-
Index: trunk/mars/libf/phymars/meso_inifis.F
===================================================================
--- trunk/mars/libf/phymars/meso_inifis.F	(revision 85)
+++ trunk/mars/libf/phymars/meso_inifis.F	(revision 86)
@@ -106,4 +106,6 @@
       CHARACTER ch1*12
       CHARACTER ch80*80
+
+#ifdef MESOSCALE
 
 !      logical chem, h2o
@@ -807,3 +809,6 @@
 !
 !      RETURN
+
+#endif
+
       END
Index: trunk/mars/libf/phymars/meso_newcondens.F
===================================================================
--- trunk/mars/libf/phymars/meso_newcondens.F	(revision 85)
+++ 	(revision )
@@ -1,877 +1,0 @@
-      SUBROUTINE meso_newcondens(ngrid,nlayer,nq,ptimestep,
-     $                  pcapcal,pplay,pplev,ptsrf,pt,
-     $                  pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,
-     $                  piceco2,psolaralb,pemisurf,
-     $                  pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc,
-     $                  fluxsurf_sw,zls)
-                                                   
-       IMPLICIT NONE
-c=======================================================================
-c   subject:
-c   --------
-c   Condensation/sublimation of CO2 ice on the ground and in the
-c   atmosphere 
-c   (Scheme described in Forget et al., Icarus, 1998)
-c
-c   author:   Francois Forget     1994-1996
-c   ------
-c
-c   input:
-c   ----- 
-c    ngrid                 nombre de points de verticales
-c                          (toutes les boucles de la physique sont au
-c                          moins vectorisees sur ngrid)
-c    nlayer                nombre de couches
-c    pplay(ngrid,nlayer)   Pressure levels 
-c    pplev(ngrid,nlayer+1) Pressure levels 
-c    pt(ngrid,nlayer)      temperature (en K)
-c    ptsrf(ngrid)          temperature de surface
-c 
-c                    \
-c    pdt(ngrid,nlayermx)   \  derivee temporelle physique avant condensation
-c                     /  ou sublimation pour  pt,ptsrf
-c    pdtsrf(ngrid)   /
-c
-c   output:
-c   -------
-c
-c    pdpsrf(ngrid)   \  derivee temporelle physique (contribution de
-c    pdtc(ngrid,nlayermx) /  la condensation ou sublimation) pour Ps,pt,ptsrf
-c    pdtsrfc(ngrid) /
-c
-c   Entree/sortie
-c   -------------
-c    
-c    piceco2(ngrid) :      quantite de glace co2 au sol (kg/m2)
-c    psolaralb(ngrid,2) :  albedo au sol
-c    pemisurf(ngrid)     :  emissivite du sol             
-
-c
-c=======================================================================
-c
-c    0.  Declarations :
-c    ------------------
-c
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comcstfi.h"
-#include "surfdat.h"
-#include "comgeomfi.h"
-#include "comvert.h"
-#include "paramet.h"
-#include "callkeys.h"
-#include "tracer.h"
-
-c-----------------------------------------------------------------------
-c    Arguments :
-c    ---------
-      INTEGER ngrid, nlayer, nq
-
-      REAL ptimestep 
-      REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
-      REAL pcapcal(ngrid)
-      REAL pt(ngrid,nlayer)
-      REAL ptsrf(ngrid)
-      REAL pphi(ngrid,nlayer)
-      REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer)
-      REAL pdtsrfc(ngrid),pdpsrf(ngrid)
-      REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid)
-
-      REAL pu(ngrid,nlayer) ,  pv(ngrid,nlayer)
-      REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
-      REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
-      REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
-      REAL pdqc(ngrid,nlayer,nq)
-      REAL fluxsurf_sw(ngrid,2) ! added to calculate flux dependent albedo 
-      real zls ! solar longitude (rad)
-c
-c    Local variables :
-c    -----------------
-
-c   variables used for albedo parametrization       
-c --------------------------------------------      
-      INTEGER i,j
-      REAL Fluxmean(jjp1)
-      INTEGER l,ig,iq,icap,nmix
-      LOGICAL transparency, fluxdependent
-c   flag transparency if you want to make the co2ice semi-transparent       
-      PARAMETER(transparency=.true.)
-c   flag fluxdependent if you want the co2ice albedo to be dependent on
-c   the incident solar flux      
-      PARAMETER(fluxdependent=.false.)
-      REAL slopy,alpha,constA,constB,constT,albediceF_new(ngridmx)
-      REAL zt(ngridmx,nlayermx)
-      REAL zcpi
-      REAL ztcond (ngridmx,nlayermx)
-      REAL ztcondsol(ngridmx) 
-      REAL zdiceco2(ngridmx)
-      REAL zcondicea(ngridmx,nlayermx)
-      REAL zcondices(ngridmx)
-      REAL zfallice(ngridmx,nlayermx+1) , zfallheat
-      REAL zmflux(nlayermx+1)
-      REAL zu(nlayermx),zv(nlayermx)
-      REAL zq(nlayermx,nqmx),zq1(nlayermx)
-      REAL ztsrf(ngridmx)
-      REAL ztc(nlayermx), ztm(nlayermx+1) 
-      REAL zum(nlayermx+1) , zvm(nlayermx+1)
-      REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
-      REAL masse(nlayermx),w(nlayermx+1)
-      REAL Sm(nlayermx),Smq(nlayermx,nqmx),mixmas,qmix
-      LOGICAL condsub(ngridmx)
-
-c variable speciale diagnostique
-      real tconda1(ngridmx,nlayermx)
-      real tconda2(ngridmx,nlayermx)
-c     REAL zdiceco2a(ngridmx) ! for diagnostic only
-      real zdtsig (ngridmx,nlayermx)
-      real zdt (ngridmx,nlayermx)
-      real vmr_co2(ngridmx,nlayermx) ! co2 volume mixing ratio
-! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer)
-! then condensation temperature is computed using partial pressure of CO2
-      logical,parameter :: improved_ztcond=.true.
-! Bound co2 (tracer) values...
-      logical,parameter :: bound_qco2=.false.
-      real,parameter :: qco2max=1.1
-      real,parameter :: qco2mini=0.1
-      real :: zqco2
-
-c   local saved variables
-      integer ico2
-      real qco2min,qco2,mmean
-      save ico2,qco2min
-      REAL emisref(ngridmx)
-      REAL latcond,tcond1mb
-      REAL acond,bcond,ccond,cpice
-!      REAL albediceF(ngridmx)
-      SAVE emisref
-      SAVE latcond,acond,bcond,ccond,cpice
-!      SAVE albediceF
-      real m_co2, m_noco2, A , B
-      save A, B, m_co2, m_noco2
-
-      LOGICAL firstcall !,firstcall2
-      SAVE firstcall !,firstcall2
-!      REAL SSUM
-!      EXTERNAL SSUM
-
-      DATA latcond,tcond1mb/5.9e5,136.27/
-      DATA cpice /1000./
-      DATA firstcall/.true./
-!      DATA firstcall2/.true./
-
-      integer flag
-
-c----------------------------------------------------------------------
-
-c   Initialisation
-c   --------------
-c
-      IF (firstcall) THEN
-         bcond=1./tcond1mb
-         ccond=cpp/(g*latcond)
-         acond=r/latcond
-
-         firstcall=.false.
-         write(*,*) 'Newcondens: improved_ztcond=',improved_ztcond
-         write(*,*) 'Newcondens: bound_qco2=',bound_qco2
-         PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
-         PRINT*,'acond,bcond,ccond',acond,bcond,ccond
-
-         ico2=0
-
-         if (tracer) then
-c          Prepare Special treatment if one of the tracer is CO2 gas
-           do iq=1,nqmx
-             if (noms(iq).eq."co2") then
-                ico2=iq
-                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
-                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
-c               Compute A and B coefficient use to compute
-c               mean molecular mass Mair defined by
-c               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
-c               1/Mair = A*q(ico2) + B
-                A =(1/m_co2 - 1/m_noco2)
-                B=1/m_noco2
-             endif
-           enddo
-c          minimum CO2 mix. ratio below which mixing occur with layer above: 
-           qco2min =0.75  
-         end if
-      ENDIF
-      zcpi=1./cpp
-c
-c======================================================================
-c    Calcul of CO2 condensation sublimation 
-c    ============================================================
-c  
-c    Used variable :
-c       piceco2(ngrid)   :  amount of co2 ice on the ground (kg/m2)
-c       zcondicea(ngrid,l):  condensation rate in layer  l  (kg/m2/s)
-c       zcondices(ngrid):  condensation rate on the ground (kg/m2/s)
-c       zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s)
-c                           
-c       pdtc(ngrid,nlayermx)  : dT/dt due to cond/sub
-c
-c
-c     Tendencies set to 0 (except pdtc)
-c     -------------------------------------
-      DO l=1,nlayer
-         DO ig=1,ngrid
-           zcondicea(ig,l) = 0.
-           zfallice(ig,l) = 0.
-           pduc(ig,l)  = 0
-           pdvc(ig,l)  = 0
-         END DO
-      END DO
-         
-      DO iq=1,nqmx         
-      DO l=1,nlayer
-         DO ig=1,ngrid
-           pdqc(ig,l,iq)  = 0
-         END DO
-      END DO
-      END DO
-
-      DO ig=1,ngrid
-         zfallice(ig,nlayer+1) = 0.
-         zcondices(ig) = 0.
-         pdtsrfc(ig) = 0.
-         pdpsrf(ig) = 0.
-         condsub(ig) = .false.
-         zdiceco2(ig) = 0.
-      ENDDO
-      zfallheat=0
-
-c     *************************
-c     ATMOSPHERIC CONDENSATION
-c     *************************
-
-c     Compute CO2 Volume mixing ratio
-c     -------------------------------
-      if (improved_ztcond.and.(ico2.ne.0)) then
-         DO l=1,nlayer
-            DO ig=1,ngrid
-              qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
-c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
-              mmean=1/(A*qco2 +B)
-              vmr_co2(ig,l) = qco2*mmean/m_co2 
-            ENDDO
-         ENDDO
-      else
-         DO l=1,nlayer
-            DO ig=1,ngrid
-              vmr_co2(ig,l)=0.95
-            ENDDO
-         ENDDO
-      end if
-
-c     forecast of atmospheric temperature zt and frost temperature ztcond
-c     --------------------------------------------------------------------
-
-      DO l=1,nlayer
-         DO ig=1,ngrid
-            zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
-!            ztcond(ig,l)=1./(bcond-acond*log(.0095*pplay(ig,l)))
-            ztcond(ig,l)=
-     &         1./(bcond-acond*log(.01*vmr_co2(ig,l)*pplay(ig,l)))
-            if (pplay(ig,l).lt.1e-4) ztcond(ig,l)=0.0 !mars Monica
-         ENDDO
-      ENDDO
- 
-c      Condensation/sublimation in the atmosphere
-c      ------------------------------------------
-c      (calcul of zcondicea , zfallice and pdtc)
-c
-      DO l=nlayer , 1, -1
-         DO ig=1,ngrid
-           pdtc(ig,l)=0.
-           IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
-               condsub(ig)=.true.
-               IF (zfallice(ig,l+1).gt.0) then  
-                 zfallheat=zfallice(ig,l+1)*
-     &           (pphi(ig,l+1)-pphi(ig,l) +
-     &          cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond
-               ELSE
-                   zfallheat=0.
-               ENDIF
-               pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
-               zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
-     &                        *ccond*pdtc(ig,l)- zfallheat
-c              Case when the ice from above sublimes entirely
-c              """""""""""""""""""""""""""""""""""""""""""""""
-               IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then
-                  pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/
-     &              (ccond*(pplev(ig,l)-pplev(ig,l+1)))
-                  zcondicea(ig,l)= -zfallice(ig,l+1)
-               END IF
-
-               zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
-            END IF
-         ENDDO
-      ENDDO
-
-c     *************************
-c       SURFACE  CONDENSATION
-c     *************************
-
-c     forecast of ground temperature ztsrf and frost temperature ztcondsol
-c     --------------------------------------------------------------------
-      DO ig=1,ngrid
-         ztcondsol(ig)=
-     &          1./(bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1)))
-         ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
-      ENDDO
-
-c
-c      Condensation/sublimation on the ground 
-c      --------------------------------------
-c      (calcul of zcondices , pdtsrfc)
-c
-      DO ig=1,ngrid
-         IF(ig.GT.ngrid/2+1) THEN
-            icap=2
-         ELSE
-            icap=1
-         ENDIF
-       
-c
-c        Loop on where we have  condensation/ sublimation 
-         IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.   ! ground cond 
-     $       (zfallice(ig,1).NE.0.) .OR.           ! falling snow
-     $      ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  ! ground sublim.
-     $      ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
-            condsub(ig) = .true. 
-
-            IF (zfallice(ig,1).gt.0) then  
-                 zfallheat=zfallice(ig,1)*
-     &           (pphi(ig,1)- phisfi(ig) +
-     &           cpice*(ztcond(ig,1)-ztcondsol(ig)))/(latcond*ptimestep)
-            ELSE
-                 zfallheat=0.
-            ENDIF
-
-c           condensation or partial sublimation of CO2 ice
-c           """""""""""""""""""""""""""""""""""""""""""""""
-            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) 
-     &      /(latcond*ptimestep)         - zfallheat
-            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
-       
-c           If the entire CO_2 ice layer sublimes
-c           """"""""""""""""""""""""""""""""""""""""""""""""""""
-c           (including what has just condensed in the atmosphere)
-
-            IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
-     &          -zcondices(ig))THEN
-               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) 
-               pdtsrfc(ig)=(latcond/pcapcal(ig))*
-     &         (zcondices(ig)+zfallheat)
-            END IF
-
-c           Changing CO2 ice amount and pressure :
-c           """"""""""""""""""""""""""""""""""""
-
-            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
-            piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
-            pdpsrf(ig) = -zdiceco2(ig)*g
-
-            IF(ABS(pdpsrf(ig)*ptimestep).GT.pplev(ig,1)) THEN
-               PRINT*,'STOP in condens'
-               PRINT*,'condensing more than total mass'
-               PRINT*,'Grid point ',ig
-               PRINT*,'Ps = ',pplev(ig,1)
-               PRINT*,'d Ps = ',pdpsrf(ig)
-               STOP
-            ENDIF
-         END IF ! if there is condensation/sublimmation
-      ENDDO ! of DO ig=1,ngrid
-
-c  ********************************************************************
-c  Surface albedo and emissivity of the surface below the snow (emisref)
-c  ********************************************************************
-c     Prepare the case where albedo varies with insolation:
-c     ----------------------------------------------------
-!      if (fluxdependent) then
-!
-c       Calcul du flux moyen (zonal mean)
-!        do j=1,jjp1
-!           Fluxmean(j)=0      
-!           do i=1,iim
-!               ig=1+(j-2)*iim +i
-!               if(j.eq.1) ig=1
-!               if(j.eq.jjp1) ig=ngrid
-!               Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1)
-!     $            +fluxsurf_sw(ig,2)
-!           enddo
-!           Fluxmean(j)=Fluxmean(j)/float(iim)
-!        enddo
-!
-c       const A and B used to calculate the albedo which depends on solar flux 
-c       albedice=constA+constB*Flux
-c       constT = time step to calculate the solar flux when flux decreases         
-!         constA=0.26
-c        constA=0.33
-c        constA=0.186
-!         constB=0.00187
-!         constT=10
-!      endif ! of if (fluxdependent)
-
-!     Check that amont of CO2 ice is not problematic
-      DO ig=1,ngrid
-           if(.not.piceco2(ig).ge.0.) THEN
-              if(piceco2(ig).le.-5.e-8) print*,
-     $         'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
-               piceco2(ig)=0.
-           endif
-      ENDDO
-      
-!     Set albedo and emissivity of the surface
-!     ----------------------------------------
-      CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
-
-c     Calcul de l'albedo
-c     ------------------
-!         do ig =1,ngrid
-!           IF(ig.GT.ngrid/2+1) THEN
-!              icap=2
-!           ELSE
-!              icap=1
-!           ENDIF
-!           IF(firstcall2) THEN
-!           albediceF(ig)=albedice(icap)
-!           ENDIF
-c       if there is still co2ice        ccccccccccccccccccccccc
-!           if (piceco2(ig).gt.0) then
-!              emisref(ig) = emisice(icap)
-
-c     if flux dependent albedo is used 
-c     --------------------------------
-!             if (fluxdependent) then
-!                j=INT((ig-2)/iim)+2
-!                if(ig.eq.1) j=1
-!                if(ig.eq.ngrid) j=jjp1
-c                albediceF_new(ig)=MIN(constA+constB*Fluxmean(j),
-c     $          constA+constB*250)
-!                  albediceF_new(ig)=constA+constB*Fluxmean(j)
-!                if (albediceF(ig).gt.albediceF_new(ig)) then
-!                   albediceF(ig)=albediceF(ig)+ ptimestep/(daysec*
-!     $             constT)*(albediceF_new(ig)-albediceF(ig))
-!                else
-!                   albediceF(ig)=albediceF_new(ig)
-!                endif
-c               if part of the ice is transparent
-c               slopy = pente de la droite:  alpha = y*co2ice/1620
-c               pour une valeur superieur a une epaisseur de glace donnee
-c               ici, epaisseur limite = 10cm
-!                if (transparency) then
-!                slopy=1/(1620*0.10)
-!                alpha=MIN(slopy*piceco2(ig),1.)
-!                psolaralb(ig,1) = alpha*albediceF(ig)
-!     $                           +(1-alpha)*albedodat(ig) 
-!                psolaralb(ig,2) = psolaralb(ig,1) 
-!                else
-!                psolaralb(ig,1) = albediceF(ig) 
-!                psolaralb(ig,2) = psolaralb(ig,1) 
-!                endif
-!             else
-c           transparency set to true and fluxdependent set to false            
-!                if (transparency) then
-!                slopy=1/(1620*0.10)
-!                alpha=MIN(slopy*piceco2(ig),1.)
-!                psolaralb(ig,1) = alpha*albedice(icap)
-!     $                           +(1-alpha)*albedodat(ig) 
-!                psolaralb(ig,2) = psolaralb(ig,1) 
-!                else
-c           simplest case: transparency and flux dependent set to false
-!                psolaralb(ig,1) = albedice(icap)
-!                psolaralb(ig,2) = albedice(icap)
-!                endif
-!             endif    
-c         no more co2ice, albedo = ground albedo
-!           else
-!              psolaralb(ig,1) = albedodat(ig)
-!              psolaralb(ig,2) = albedodat(ig)
-!              emisref(ig) = emissiv
-!              pemisurf(ig) = emissiv
-!           endif
-!         end do ! end of the ig loop
-
-! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
-      DO ig=1,ngrid
-        if (piceco2(ig).eq.0) then
-          pemisurf(ig)=emissiv
-        endif
-      ENDDO
-
-!         firstcall2=.false.
-c ***************************************************************
-c  Correction to account for redistribution between sigma or hybrid 
-c  layers when changing surface pressure (and warming/cooling
-c  of the CO2 currently changing phase).
-c *************************************************************
-
-      DO ig=1,ngrid
-        if (condsub(ig)) then
-           do l=1,nlayer
-             ztc(l)  =zt(ig,l)   +pdtc(ig,l)  *ptimestep
-             zu(l)   =pu(ig,l)   +pdu( ig,l)  *ptimestep
-             zv(l)   =pv(ig,l)   +pdv( ig,l)  *ptimestep
-            do iq=1,nqmx
-             zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
-            enddo
-           end do
-
-c  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
-c  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
-
-            zmflux(1) = -zcondices(ig)
-            DO l=1,nlayer
-             zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
-     &        + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
-c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
-          if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
-            END DO
-
-c Mass of each layer
-c ------------------ 
-            DO l=1,nlayer
-              masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
-            END DO
-
-
-c  Corresponding fluxes for T,U,V,Q
-c  """"""""""""""""""""""""""""""""
-
-c           averaging operator for TRANSPORT  
-c           """"""""""""""""""""""""""""""""
-c           Value transfert at the surface interface when condensation
-c           sublimation:
-            ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
-            zum(1) = 0  
-            zvm(1) = 0  
-            do iq=1,nqmx
-              zqm(1,iq)=0. ! most tracer do not condense !
-            enddo
-c           Special case if one of the tracer is CO2 gas
-            if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2
-
-c           Van Leer scheme:
-            DO l=1,nlayer+1
-                w(l)=-zmflux(l)*ptimestep
-            END DO
-            call vl1d(ztc,2.,masse,w,ztm) 
-            call vl1d(zu ,2.,masse,w,zum) 
-            call vl1d(zv ,2.,masse,w,zvm) 
-            do iq=1,nqmx
-             do l=1,nlayer
-              zq1(l)=zq(l,iq)
-             enddo
-             zqm1(1)=zqm(1,iq)
-             call vl1d(zq1,2.,masse,w,zqm1)
-             do l=2,nlayer
-              zq( l,iq)=zq1(l)
-              zqm(l,iq)=zqm1(l)
-             enddo
-            enddo
-
-c           Surface condensation affects low winds
-            if (zmflux(1).lt.0) then 
-                zum(1)= zu(1) *  (w(1)/masse(1))
-                zvm(1)= zv(1) *  (w(1)/masse(1))
-                if (w(1).gt.masse(1)) then ! ensure numerical stability
-                  zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2)
-                  zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2)
-                end if
-            end if
-                    
-            ztm(nlayer+1)= ztc(nlayer) ! should not be used, but... 
-            zum(nlayer+1)= zu(nlayer)  ! should not be used, but...
-            zvm(nlayer+1)= zv(nlayer)  ! should not be used, but...
-            do iq=1,nqmx
-             zqm(nlayer+1,iq)= zq(nlayer,iq)
-            enddo
- 
-CCCC
-CCCC *** WRF comments
-CCCC
-cc           Tendencies on T, U, V, Q 
-cc           """"""""""""""""""""""""
-!            DO l=1,nlayer
-! 
-!c             Tendencies on T
-!                zdtsig(ig,l) = (1/masse(l)) *
-!     &        ( zmflux(l)*(ztm(l) - ztc(l)) 
-!     &        - zmflux(l+1)*(ztm(l+1) - ztc(l))
-!     &        + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )
-!                pdtc(ig,l) =  pdtc(ig,l) + zdtsig(ig,l)
-!
-!c             Tendencies on U
-!                pduc(ig,l)   = (1/masse(l)) *
-!     &        ( zmflux(l)*(zum(l) - zu(l))
-!     &        - zmflux(l+1)*(zum(l+1) - zu(l)) )
-!
-!
-!c             Tendencies on V
-!                pdvc(ig,l)   = (1/masse(l)) *
-!     &        ( zmflux(l)*(zvm(l) - zv(l))
-!     &        - zmflux(l+1)*(zvm(l+1) - zv(l)) )
-!
-!            END DO
-
-c           Tendencies on Q
-            do iq=1,nqmx
-!              if (noms(iq).eq.'co2') then 
-              if (iq.eq.ico2) then 
-c               SPECIAL Case when the tracer IS CO2 :
-                DO l=1,nlayer
-                  pdqc(ig,l,iq)= (1/masse(l)) *
-     &           ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
-     &           - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) 
-     &           + zcondicea(ig,l)*(zq(l,iq)-1.) )
-                END DO
-              else
-                DO l=1,nlayer
-                  pdqc(ig,l,iq)= (1/masse(l)) *
-     &           ( zmflux(l)*(zqm(l,iq) - zq(l,iq))
-     &           - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) 
-     &           + zcondicea(ig,l)*zq(l,iq) )
-                END DO
-              end if
-            enddo
-
-c           --------------------------------------------------------
-c           Roughly Simulate Molecular mixing when CO2 is too depleted by
-c           Surface condensation (mixing starts if qco2 < qco2min )  
-c           FF 06/2004
-c           WARNING : this is now done in convadj, better (FF 02/2005)
-c           --------------------------------------------------------
-            flag=0  ! now done in convadj : must be =0
-            if (flag.eq.1) then 
-            if(ico2.gt.0) then   ! relevant only if one tracer is CO2
-             if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep
-     &       .lt.qco2min) then 
-                do iq=1,nqmx
-                  zq(1,iq)=pq(ig,1,iq)
-     &                     +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep
-                  Smq(1,iq) = masse(1)*zq(1,iq)
-                end do
-                Sm(1)  = masse(1)
-                do l =2,nlayermx
-                  do iq=1,nqmx
-                     zq(l,iq)=pq(ig,l,iq)
-     &                        +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep
-                     smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq)
-                  end do 
-                  sm(l) = sm(l-1) + masse(l)
-                  if(zq(l,ico2).gt.qco2min) then 
-c                   mixmas: mass of atmosphere that must be mixed to reach qco2min
-                    mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2))
-     &                      /(zq(l,ico2)-qco2min)
-                    if((mixmas.le.sm(l)))then
-c                      OK if mixed mass less than mass of layers affected
-                       nmix=l   ! number of layer affected by mixing
-                       goto 99 
-                    end if
-                  end if
-                 end do
- 99              continue 
-                 do iq=1,nqmx
-                   qmix=zq(nmix,iq) 
-     &             +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas
-                   do l=1,nmix-1
-                      pdqc(ig,l,iq)=
-     &                  (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq)
-                   end do 
-c                  layer only partly mixed :
-                   pdqc(ig,nmix,iq)=(
-     &             qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix)
-     &             -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq)
-
-                 end do 
-             end if
-            end if
-
-        endif   ! (flag.eq.1)
-       end if ! if (condsub)
-      END DO  ! loop on ig 
-
-c ***************************************************************
-c   CO2 snow / clouds scheme
-c ***************************************************************
-
-      call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
-     &        zcondicea,zcondices,zfallice,pemisurf)
-
-c ***************************************************************
-c Ecriture des diagnostiques
-c ***************************************************************
-
-c     DO l=1,nlayer
-c        DO ig=1,ngrid
-c         Taux de cond en kg.m-2.pa-1.s-1
-c          tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
-c          Taux de cond en kg.m-3.s-1
-c          tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l))
-c        END DO
-c     END DO
-c     call WRITEDIAGFI(ngridmx,'tconda1',
-c    &'Taux de condensation CO2 atmospherique /Pa',
-c    & 'kg.m-2.Pa-1.s-1',3,tconda1)
-c     call WRITEDIAGFI(ngridmx,'tconda2',
-c    &'Taux de condensation CO2 atmospherique /m',
-c    & 'kg.m-3.s-1',3,tconda2)
-
-! output falling co2 ice in 1st layer:
-!      call WRITEDIAGFI(ngridmx,'fallice',
-!     &'Precipitation of co2 ice',
-!     & 'kg.m-2.s-1',2,zfallice(1,1))
-
-!! Specific stuff to bound co2 tracer values ....
-      if (bound_qco2.and.(ico2.ne.0)) then
-        do ig=1,ngridmx
-          do l=1,nlayermx
-            zqco2=pq(ig,l,ico2)
-     &            +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep
-            if (zqco2.gt.qco2max) then
-            ! correct pdqc:
-              pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep)
-     &                               -pdq(ig,l,ico2)
-              write(*,*) "newcondens: adapting pdqc(ig,l,ico2)",
-     &                   " so that co2 conc. does not exceed",qco2max
-              write(*,*) "   ig:",ig,"  l:",l
-            endif ! of if (zqco2.gt.qco2max)
-            if (zqco2.lt.qco2mini) then
-            ! correct pdqc:
-              pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep)
-     &                               -pdq(ig,l,ico2)
-              write(*,*) "newcondens: adapting pdqc(ig,l,ico2)",
-     &                   " so that co2 conc. is not less than",qco2mini
-              write(*,*) "   ig:",ig,"  l:",l
-            endif ! of if (zqco2.lt.qco2mini)
-          end do 
-        enddo
-      endif ! of if (bound_qco2.and.(ico2.ne.0)) then
-
-      return
-      end
-
-
-
-c *****************************************************************
-      SUBROUTINE vl1d(q,pente_max,masse,w,qm)
-c
-c    
-c     Operateur de moyenne inter-couche pour calcul de transport type
-c     Van-Leer " pseudo amont " dans la verticale
-c    q,w sont des arguments d'entree  pour le s-pg ....
-c    masse : masse de la couche Dp/g
-c    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
-c    pente_max = 2 conseillee
-c
-c
-c   --------------------------------------------------------------------
-      IMPLICIT NONE
-
-#include "dimensions.h"
-
-c
-c
-c
-c   Arguments:
-c   ----------
-      real masse(llm),pente_max
-      REAL q(llm),qm(llm+1)
-      REAL w(llm+1)
-c
-c      Local 
-c   ---------
-c
-      INTEGER l
-c
-      real dzq(llm),dzqw(llm),adzqw(llm),dzqmax
-      real sigw, Mtot, MQtot
-      integer m 
-c     integer ismax,ismin 
-
-
-c    On oriente tout dans le sens de la pression 
-c     W > 0 WHEN DOWN !!!!!!!!!!!!!
-
-      do l=2,llm
-            dzqw(l)=q(l-1)-q(l)
-            adzqw(l)=abs(dzqw(l))
-      enddo
-
-      do l=2,llm-1
-            if(dzqw(l)*dzqw(l+1).gt.0.) then
-                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
-            else
-                dzq(l)=0.
-            endif
-            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
-            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
-      enddo
-
-         dzq(1)=0.
-         dzq(llm)=0.
-
-       do l = 1,llm-1
-
-c         Regular scheme (transfered mass < layer mass)
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
-             sigw=w(l+1)/masse(l+1)
-             qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
-          else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
-             sigw=w(l+1)/masse(l)
-             qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))
-
-c         Extended scheme (transfered mass > layer mass)
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          else if(w(l+1).gt.0.) then
-             m=l+1
-             Mtot = masse(m)
-             MQtot = masse(m)*q(m)
-             do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1))))
-                m=m+1
-                Mtot = Mtot + masse(m)
-                MQtot = MQtot + masse(m)*q(m)
-             end do
-             if (m.lt.llm) then
-                sigw=(w(l+1)-Mtot)/masse(m+1)
-                qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)*
-     &          (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
-             else
-                w(l+1) = Mtot
-                qm(l+1) = Mqtot / Mtot
-                write(*,*) 'top layer is disapearing !'
-                stop
-             end if
-          else      ! if(w(l+1).lt.0) 
-             m = l-1 
-             Mtot = masse(m+1)
-             MQtot = masse(m+1)*q(m+1)
-             do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
-                m=m-1
-                Mtot = Mtot + masse(m+1)
-                MQtot = MQtot + masse(m+1)*q(m+1)
-             end do
-             if (m.gt.0) then
-                sigw=(w(l+1)+Mtot)/masse(m)
-                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* 
-     &          (q(m)-0.5*(1.+sigw)*dzq(m))  )
-             else
-                qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
-             end if   
-          end if
-       enddo
-
-c     boundary conditions (not used in newcondens !!)
-c         qm(llm+1)=0.
-c         if(w(1).gt.0.) then
-c            qm(1)=q(1)
-c         else 
-c           qm(1)=0.
-c         end if
-
-      return
-      end
Index: trunk/mars/libf/phymars/meso_physiq.F
===================================================================
--- trunk/mars/libf/phymars/meso_physiq.F	(revision 85)
+++ trunk/mars/libf/phymars/meso_physiq.F	(revision 86)
@@ -10,5 +10,7 @@
      $            wisoil,wdsoil,
      $            wecritphys,
+#ifdef MESOSCALE
      $            output_tab2d, output_tab3d,
+#endif
      $            flag_LES)
 
@@ -152,7 +154,10 @@
 #include "netcdf.inc"
 
+!!!!**** SPECIFIC TO MESOSCALE
+#ifdef MESOSCALE
 #include "meso_slope.h"
 #include "wrf_output_2d.h"
 #include "wrf_output_3d.h"
+#endif
 
 #include "advtrac.h"   !!! this is necessary for tracers (in dyn3d)
@@ -172,5 +177,5 @@
       REAL zh(ngridmx,nlayermx)      ! potential temperature (K)
       LOGICAL firstcall,lastcall
-c ****WRF
+!!! ****WRF WRF specific to mesoscale
       INTEGER wday_ini
       REAL wtsurf(ngridmx)  ! input only ay firstcall - output
@@ -183,6 +188,8 @@
       REAL wq2(ngridmx,nlayermx+1)
       REAL wecritphys
+#ifdef MESOSCALE
       REAL output_tab2d(ngridmx,n2d)
       REAL output_tab3d(ngridmx,nlayer,n3d)
+#endif
       REAL sl_ls, sl_lct, sl_lat, sl_tau, sl_alb, sl_the, sl_psi
       REAL sl_fl0, sl_flu
@@ -196,5 +203,5 @@
       integer iloop
       INTEGER tracerset    !!! this corresponds to config%mars
-c ****WRF
+!!! ****WRF WRF specific to mesoscale
       REAL pday
       REAL ptime 
@@ -353,4 +360,5 @@
 
 c=======================================================================
+#ifdef MESOSCALE
 
 c 1. Initialisation:
@@ -1036,5 +1044,5 @@
 
       IF (callcond) THEN
-         CALL meso_newcondens(ngrid,nlayer,nq,ptimestep,
+         CALL newcondens(ngrid,nlayer,nq,ptimestep,
      $              capcal,pplay,pplev,tsurf,pt,
      $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
@@ -1890,4 +1898,5 @@
       icount=icount+1
       write(*,*) 'now, back to the dynamical core...'
+#endif
       RETURN
       END
Index: trunk/mars/libf/phymars/meso_testphys1d.F
===================================================================
--- trunk/mars/libf/phymars/meso_testphys1d.F	(revision 85)
+++ 	(revision )
@@ -1,687 +1,0 @@
-
-      PROGRAM meso_testphys1d
-! to use  'getin'
-      USE ioipsl_getincom 
-      IMPLICIT NONE
-
-c=======================================================================
-c   subject:
-c   --------
-c   PROGRAM useful to run physical part of the martian GCM in a 1D column
-c       
-c Can be compiled with a command like (e.g. for 25 layers)
-c  "makegcm -p mars -d 25 testphys1d"
-c It requires the files "testphys1d.def" "callphys.def"
-c   and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line)
-c      and a file describing the sigma layers (e.g. "z2sig.def")
-c
-c   author: Frederic Hourdin, R.Fournier,F.Forget
-c   -------
-c   
-c   update: 12/06/2003 including chemistry (S. Lebonnois) 
-c                            and water ice (F. Montmessin)
-c 
-c=======================================================================
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "dimradmars.h"
-#include "comgeomfi.h"
-#include "surfdat.h"
-#include "comsoil.h"
-#include "comdiurn.h"
-#include "callkeys.h"
-#include "comcstfi.h"
-#include "planete.h"
-#include "comsaison.h"
-#include "yomaer.h"
-#include "control.h"
-#include "comvert.h"
-#include "netcdf.inc"
-#include "comg1d.h"
-#include "watercap.h"
-#include "logic.h"
-#include "advtrac.h"
-
-c --------------------------------------------------------------
-c  Declarations
-c --------------------------------------------------------------
-c
-      INTEGER unitstart      ! unite d'ecriture de "startfi"
-      INTEGER nlayer,nlevel,nsoil,ndt
-      INTEGER ilayer,ilevel,isoil,idt,iq
-      LOGICAl firstcall,lastcall
-c
-      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
-      REAL day           ! date durant le run
-      REAL time             ! time (0<time<1 ; time=0.5 a midi)
-      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
-      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
-      REAL psurf,tsurf      
-      REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
-      REAL gru,grv   ! prescribed "geostrophic" background wind
-      REAL temp(nlayermx)   ! temperature at the middle of the layers
-      REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
-      REAL qsurf(nqmx)      ! tracer surface budget (e.g. kg.m-2)
-      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
-      REAL co2ice           ! co2ice layer (kg.m-2)
-      REAL emis             ! surface layer
-      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
-      REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
-
-c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
-      REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
-      REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
-      REAL dpsurf    
-      REAL dq(nlayermx,nqmx)
-      REAL dqdyn(nlayermx,nqmx)
-
-c   Various intermediate variables
-      INTEGER thermo
-      REAL zls
-      REAL phi(nlayermx),h(nlayermx),s(nlayermx)
-      REAL pks, ptif, w(nlayermx)
-      REAL qtotinit, mqtot(nqmx),qtot
-      INTEGER ierr, aslun
-      REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
-      Logical  tracerdyn
-      integer :: nq=1 ! number of tracers
-
-      character*2 str2
-      character (len=7) :: str7
-      character(len=44) :: txt
-
-c=======================================================================
-
-c=======================================================================
-c INITIALISATION
-c=======================================================================
-
-c ------------------------------------------------------
-c  Constantes prescrites ICI
-c ------------------------------------------------------
-
-      pi=2.E+0*asin(1.E+0)
-
-c     Constante de la Planete Mars
-c     ----------------------------
-      rad=3397200.               ! rayon de mars (m)  ~3397200 m
-      daysec=88775.              ! duree du sol (s)  ~88775 s
-      omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
-      g=3.72                     ! gravite (m.s-2) ~3.72  
-      mugaz=43.49                ! Masse molaire de l'atm (g.mol-1) ~43.49
-      rcp=.256793                ! = r/cp  ~0.256793
-      r= 8.314511E+0 *1000.E+0/mugaz
-      cpp= r/rcp
-      year_day = 669           ! duree de l'annee (sols) ~668.6
-      periheli = 206.66          ! dist.min. soleil-mars (Mkm) ~206.66
-      aphelie = 249.22           ! dist.max. soleil-mars (Mkm) ~249.22
-      peri_day =  485.           ! date du perihelie (sols depuis printemps)
-      obliquit = 25.2           ! Obliquite de la planete (deg) ~25.2         
- 
-c     Parametres Couche limite et Turbulence 
-c     --------------------------------------
-      z0 =  1.e-2                ! surface roughness (m) ~0.01 
-      emin_turb = 1.e-6          ! energie minimale ~1.e-8
-      lmixmin = 30               ! longueur de melange ~100
- 
-c     propriete optiques des calottes et emissivite du sol
-c     ----------------------------------------------------
-      emissiv= 0.95              ! Emissivite du sol martien ~.95
-      emisice(1)=0.95            ! Emissivite calotte nord
-      emisice(2)=0.95            ! Emissivite calotte sud  
-      albedice(1)=0.5           ! Albedo calotte nord
-      albedice(2)=0.5            ! Albedo calotte sud
-      iceradius(1) = 100.e-6     ! mean scat radius of CO2 snow (north)
-      iceradius(2) = 100.e-6     ! mean scat radius of CO2 snow (south)
-      dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
-      dtemisice(2) = 2.          ! time scale for snow metamorphism (south
-      hybrid=.false.
-
-c ------------------------------------------------------
-c  Lecture des parametres dans "run.def" 
-c ------------------------------------------------------
-
-
-! check if 'run.def' file is around (otherwise reading parameters
-! from callphys.def via getin() routine won't work.
-      open(99,file='run.def',status='old',form='formatted',
-     &     iostat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) 'Cannot find required file "run.def"'
-        write(*,*) '  (which should contain some input parameters'
-        write(*,*) '   along with the following line:'
-        write(*,*) '   INCLUDEDEF=callphys.def'
-        write(*,*) '   )'
-        write(*,*) ' ... might as well stop here ...'
-        stop
-      else
-        close(99)
-      endif
-
-! check if we are going to run with or without tracers
-      write(*,*) "Run with or without tracer transport ?"
-      tracer=.false. ! default value
-      call getin("tracer",tracer)
-      write(*,*) " tracer = ",tracer
-
-! while we're at it, check if there is a 'traceur.def' file
-! and preocess it, if necessary. Otherwise initialize tracer names
-      if (tracer) then
-      ! load tracer names from file 'traceur.def'
-        open(90,file='traceur.def',status='old',form='formatted',
-     &       iostat=ierr)
-        if (ierr.ne.0) then
-          write(*,*) 'Cannot find required file "traceur.def"'
-          write(*,*) ' If you want to run with tracers, I need it'
-          write(*,*) ' ... might as well stop here ...'
-          stop
-        else
-          write(*,*) "testphys1d: Reading file traceur.def"
-          ! read number of tracers:
-          read(90,*,iostat=ierr) nq
-          if (ierr.ne.0) then
-            write(*,*) "testphys1d: error reading number of tracers"
-            write(*,*) "   (first line of traceur.def) "
-            stop
-          else
-            ! check that the number of tracers is indeed nqmx
-            if (nq.ne.nqmx) then
-              write(*,*) "testphys1d: error, wrong number of tracers:"
-              write(*,*) "nq=",nq," whereas nqmx=",nqmx
-              stop
-            endif
-          endif
-        endif
-        ! read tracer names from file traceur.def
-        do iq=1,nqmx
-          read(90,*,iostat=ierr) tnom(iq)
-          if (ierr.ne.0) then
-            write(*,*) 'testphys1d: error reading tracer names...'
-            stop
-          endif
-        enddo
-        close(90)
-
-        ! initialize tracers here:
-        write(*,*) "testphys1d: initializing tracers"
-        q(:,:)=0 ! default, set everything to zero
-        qsurf(:)=0
-        ! "smarter" initialization of some tracers
-        ! (get values from "profile_*" files, if these are available)
-        do iq=1,nqmx
-          txt=""
-          write(txt,"(a)") tnom(iq)
-          write(*,*)"  tracer:",trim(txt)
-          ! CO2
-          if (txt.eq."co2") then
-            q(:,iq)=0.95   ! kg /kg of atmosphere
-            qsurf(iq)=0. ! kg/m2 (not used for CO2)
-            ! even better, look for a "profile_co2" input file
-            open(91,file='profile_co2',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayermx
-                read(91,*) q(ilayer,iq)
-              enddo
-            endif
-            close(91)
-          endif ! of if (txt.eq."co2")
-          ! WATER VAPOUR
-          if (txt.eq."h2o_vap") then
-            !look for a "profile_h2o_vap" input file
-            open(91,file='profile_h2o_vap',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayermx
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_h2o_vap file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."h2o_ice")
-          ! WATER ICE
-          if (txt.eq."h2o_ice") then
-            !look for a "profile_h2o_vap" input file
-            open(91,file='profile_h2o_ice',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayermx
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_h2o_ice file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."h2o_ice")
-          ! DUST
-          !if (txt(1:4).eq."dust") then
-          !  q(:,iq)=0.4    ! kg/kg of atmosphere
-          !  qsurf(iq)=100 ! kg/m2
-          !endif
-          ! DUST MMR
-          if (txt.eq."dust_mass") then
-            !look for a "profile_dust_mass" input file
-            open(91,file='profile_dust_mass',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayermx
-                read(91,*) q(ilayer,iq)
-!                write(*,*) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_dust_mass file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."dust_mass")
-          ! DUST NUMBER
-          if (txt.eq."dust_number") then
-            !look for a "profile_dust_number" input file
-            open(91,file='profile_dust_number',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayermx
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_dust_number file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."dust_number")
-        enddo ! of do iq=1,nqmx
-        ! NB: some more initializations (chemistry) is done later
-
-      else
-      ! we still need to set (dummy) tracer names for physdem1
-        nq=nqmx
-        do iq=1,nq
-          write(str7,'(a1,i2.2)')'q',iq
-          tnom(iq)=str7
-        enddo
-      endif ! of if (tracer)
-
-c  Date et heure locale du debut du run
-c  ------------------------------------
-c    Date (en sols depuis le solstice de printemps) du debut du run
-      day0 = 0 ! default value for day0
-      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
-      call getin("day0",day0)
-      day=float(day0)
-      write(*,*) " day0 = ",day0
-c  Heure de demarrage
-      time=0 ! default value for time
-      write(*,*)'Initial local time (in hours, between 0 and 24)?'
-      call getin("time",time)
-      write(*,*)" time = ",time
-      time=time/24.E+0 ! convert time (hours) to fraction of sol
-
-c  Discretisation (Definition de la grille et des pas de temps)
-c  --------------
-c
-      nlayer=nlayermx
-      nlevel=nlayer+1
-      nsoil=nsoilmx
-
-      day_step=48 ! default value for day_step
-      PRINT *,'Number of time steps per sol ?'
-      call getin("day_step",day_step)
-      write(*,*) " day_step = ",day_step
-
-      ndt=10 ! default value for ndt
-      PRINT *,'Number of sols to run ?'
-      call getin("ndt",ndt)
-      write(*,*) " ndt = ",ndt
-
-      ndt=ndt*day_step     
-      dtphys=daysec/day_step  
-c Pression de surface sur la planete
-c ------------------------------------
-c
-      psurf=610. ! default value for psurf
-      PRINT *,'Surface pressure (Pa) ?'
-      call getin("psurf",psurf)
-      write(*,*) " psurf = ",psurf
-c Pression de reference
-      pa=20.
-      preff=610.      
- 
-c Proprietes des poussiere aerosol
-c --------------------------------
-      tauvis=0.2 ! default value for tauvis
-      print *,'Reference dust opacity at 700 Pa ?'
-      call getin("tauvis",tauvis)
-      write(*,*) " tauvis = ",tauvis
- 
-c  latitude/longitude
-c  ------------------
-      lati(1)=0 ! default value for lati(1)
-      PRINT *,'latitude (in degrees) ?'
-      call getin("latitude",lati(1))
-      write(*,*) " latitude = ",lati(1)
-      lati(1)=lati(1)*pi/180.E+0
-      long(1)=0.E+0
-      long(1)=long(1)*pi/180.E+0
-
-c  Initialisation albedo / inertie du sol
-c  --------------------------------------
-c
-      albedodat(1)=0.2 ! default value for albedodat
-      PRINT *,'Albedo of bare ground ?'
-      call getin("albedo",albedodat(1))
-      write(*,*) " albedo = ",albedodat(1)
-
-      inertiedat(1,1)=400 ! default value for inertiedat
-      PRINT *,'Soil thermal inertia (SI) ?'
-      call getin("inertia",inertiedat(1,1))
-      write(*,*) " inertia = ",inertiedat(1,1)
-c
-c  pour le schema d'ondes de gravite
-c  ---------------------------------
-c
-      zmea(1)=0.E+0
-      zstd(1)=0.E+0
-      zsig(1)=0.E+0
-      zgam(1)=0.E+0
-      zthe(1)=0.E+0
-
-
-
-c   Initialisation speciales "physiq"
-c   ---------------------------------
-c   la surface de chaque maille est inutile en 1D --->
-      area(1)=1.E+0
-
-c   le geopotentiel au sol est inutile en 1D car tout est controle
-c   par la pression de surface --->
-      phisfi(1)=0.E+0
-
-c  "inifis" reproduit un certain nombre d'initialisations deja faites
-c  + lecture des clefs de callphys.def
-c  + calcul de la frequence d'appel au rayonnement
-c  + calcul des sinus et cosinus des longitude latitude
-
-!Mars possible matter with dtphys in input and include!!!
-      CALL meso_inifis(1,llm,day0,daysec,dtphys,
-     .            lati,long,area,rad,g,r,cpp)
-c   Initialisation pour prendre en compte les vents en 1-D
-c   ------------------------------------------------------
-      ptif=2.E+0*omeg*sinlat(1)
- 
-c    vent geostrophique
-      gru=10. ! default value for gru
-      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
-      call getin("u",gru)
-      write(*,*) " u = ",gru
-      grv=0. !default value for grv
-      PRINT *,'meridional northward component of the geostrophic',
-     &' wind (m/s) ?'
-      call getin("v",grv)
-      write(*,*) " v = ",grv
-
-c     Initialisation des vents  au premier pas de temps
-      DO ilayer=1,nlayer
-         u(ilayer)=gru
-         v(ilayer)=grv
-      ENDDO
-
-c     energie cinetique turbulente
-      DO ilevel=1,nlevel
-         q2(ilevel)=0.E+0
-      ENDDO
-
-c  glace de CO2 au sol
-c  -------------------
-      co2ice=0.E+0 ! default value for co2ice
-      PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
-      call getin("co2ice",co2ice)
-      write(*,*) " co2ice = ",co2ice
-
-c
-c  emissivite
-c  ----------
-      emis=emissiv
-      IF (co2ice.eq.1.E+0) THEN
-         emis=emisice(1) ! northern hemisphere
-         IF(lati(1).LT.0) emis=emisice(2) ! southern hemisphere
-      ENDIF
-
- 
-
-c  calcul des pressions et altitudes en utilisant les niveaux sigma
-c  ----------------------------------------------------------------
-
-c    Vertical Coordinates
-c    """"""""""""""""""""
-      hybrid=.true.
-      PRINT *,'Hybrid coordinates ?'
-      call getin("hybrid",hybrid)
-      write(*,*) " hybrid = ", hybrid
-
-      CALL  disvert
-
-      DO ilevel=1,nlevel
-        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
-      ENDDO
-
-      DO ilayer=1,nlayer
-        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
-      ENDDO
-
-      DO ilayer=1,nlayer
-        zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1))
-     &   /g
-      ENDDO
-
-
-c  profil de temperature au premier appel
-c  --------------------------------------
-      pks=psurf**rcp
-
-c altitude en km dans profile: on divise zlay par 1000
-      tmp1(0)=0.E+0
-      DO ilayer=1,nlayer
-        tmp1(ilayer)=zlay(ilayer)/1000.E+0
-      ENDDO
-
-      call profile(nlayer+1,tmp1,tmp2)
-
-      tsurf=tmp2(0)
-      DO ilayer=1,nlayer
-        temp(ilayer)=tmp2(ilayer)
-      ENDDO
-      
-
-
-! Initialize soil properties and temperature
-! ------------------------------------------
-      volcapa=1.e6 ! volumetric heat capacity
-      DO isoil=1,nsoil
-         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
-         tsoil(isoil)=tsurf  ! soil temperature
-      ENDDO
-
-! Initialize depths
-! -----------------
-      do isoil=0,nsoil-1
-        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
-      enddo
-      do isoil=1,nsoil
-        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
-      enddo
-
-c    Initialisation des traceurs
-c    ---------------------------
-
-      if (photochem.or.callthermos) then
-         write(*,*) 'Especes chimiques initialisees'
-         ! thermo=0: initialisation pour toutes les couches 
-         thermo=0
-         call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
-     $   thermo,qsurf)
-      endif
-      watercaptag(ngridmx)=.false.
-      
-
-c    Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
-c    ----------------------------------------------------------------
-c    (effectuee avec "writeg1d" a partir de "physiq.F"
-c    ou tout autre subroutine physique, y compris celle ci).
-
-        g1d_nlayer=nlayer
-        g1d_nomfich='g1d.dat'
-        g1d_unitfich=40
-        g1d_nomctl='g1d.ctl'
-        g1d_unitctl=41
-        g1d_premier=.true.
-        g2d_premier=.true.
-
-c  Ecriture de "startfi"
-c  --------------------
-c  (Ce fichier sera aussitot relu au premier
-c   appel de "physiq", mais il est necessaire pour passer
-c   les variables purement physiques a "physiq"...
-
-      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
-     .              dtphys,float(day0),
-     .              time,tsurf,tsoil,co2ice,emis,q2,qsurf,
-     .              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
-c=======================================================================
-c  BOUCLE TEMPORELLE DU MODELE 1D 
-c=======================================================================
-c
-      firstcall=.true.
-      lastcall=.false.
-
-      DO idt=1,ndt
-c        IF (idt.eq.ndt) lastcall=.true.
-        IF (idt.eq.ndt-day_step-1) then       !test
-         lastcall=.true.
-         call solarlong(day*1.0,zls)
-         write(103,*) 'Ls=',zls*180./pi
-         write(103,*) 'Lat=', lati(1)*180./pi
-         write(103,*) 'Tau=', tauvis/700*psurf
-         write(103,*) 'RunEnd - Atmos. Temp. File'
-         write(103,*) 'RunEnd - Atmos. Temp. File'
-         write(104,*) 'Ls=',zls*180./pi
-         write(104,*) 'Lat=', lati(1)
-         write(104,*) 'Tau=', tauvis/700*psurf
-         write(104,*) 'RunEnd - Atmos. Temp. File'
-        ENDIF
-
-c    calcul du geopotentiel 
-c     ~~~~~~~~~~~~~~~~~~~~~
-      DO ilayer=1,nlayer
-        s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
-        h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
-      ENDDO
-      phi(1)=pks*h(1)*(1.E+0-s(1))
-      DO ilayer=2,nlayer
-         phi(ilayer)=phi(ilayer-1)+
-     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
-     &                  *(s(ilayer-1)-s(ilayer))
-
-      ENDDO
-
-c       appel de la physique
-c       --------------------
-
-      CALL meso_physiq (1,llm,nqmx,
-     ,     firstcall,lastcall,
-     ,     day,time,dtphys,
-     ,     plev,play,phi,
-     ,     u, v,temp, q,  
-     ,     w,
-C - sorties
-     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
-
-c       evolution du vent : modele 1D
-c       -----------------------------
- 
-c       la physique calcule les derivees temporelles de u et v.
-c       on y rajoute betement un effet Coriolis.
-c
-c       DO ilayer=1,nlayer
-c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
-c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
-c       ENDDO
-
-c       Pour certain test : pas de coriolis a l'equateur
-c       if(lati(1).eq.0.) then
-          DO ilayer=1,nlayer
-             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
-             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
-          ENDDO
-c       end if
-c      
-c
-c       Calcul du temps au pas de temps suivant
-c       ---------------------------------------
-        firstcall=.false.
-        time=time+dtphys/daysec
-        IF (time.gt.1.E+0) then
-            time=time-1.E+0
-            day=day+1
-        ENDIF
-
-c       calcul des vitesses et temperature au pas de temps suivant
-c       ----------------------------------------------------------
-
-        DO ilayer=1,nlayer
-           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
-           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
-           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
-        ENDDO
-
-c       calcul des pressions au pas de temps suivant
-c       ----------------------------------------------------------
-
-           psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
-           DO ilevel=1,nlevel
-             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
-           ENDDO
-           DO ilayer=1,nlayer
-             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
-           ENDDO
-
-!       increment tracer values
-        DO iq = 1, nqmx
-          DO ilayer=1,nlayer
-             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
-          ENDDO
-        ENDDO
-
-      ENDDO   ! fin de la boucle temporelle
-
-c    ========================================================
-c    GESTION DES SORTIE
-c    ========================================================
-
-c    fermeture pour conclure les sorties format grads dans "g1d.dat"
-c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
-c    ou tout autre subroutine physique, y compris celle ci).
-
-c        CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
-        CALL endg1d(1,nlayer,zlay/1000.,ndt)
-
-c    ========================================================
-      END
- 
-c***********************************************************************
-c***********************************************************************
-c     Subroutines Bidons utilise seulement en 3D, mais
-c     necessaire a la compilation de testphys1d en 1-D
-
-      subroutine gr_fi_dyn
-      RETURN
-      END
- 
-c***********************************************************************
-c***********************************************************************
-
-#include "../dyn3d/disvert.F"
Index: trunk/mars/libf/phymars/newcondens.F
===================================================================
--- trunk/mars/libf/phymars/newcondens.F	(revision 85)
+++ trunk/mars/libf/phymars/newcondens.F	(revision 86)
@@ -590,5 +590,10 @@
              zqm(nlayer+1,iq)= zq(nlayer,iq)
             enddo
-  
+
+#ifdef MESOSCALE
+!!!! AS: This part must be commented in the mesoscale model
+!!!! AS: ... to avoid instabilities.
+!!!! AS: you have to compile with -DMESOSCALE to do so
+#else  
 c           Tendencies on T, U, V, Q 
 c           """"""""""""""""""""""""
@@ -614,4 +619,5 @@
 
             END DO
+#endif
 
 c           Tendencies on Q
Index: trunk/mars/libf/phymars/testphys1d.F
===================================================================
--- trunk/mars/libf/phymars/testphys1d.F	(revision 85)
+++ trunk/mars/libf/phymars/testphys1d.F	(revision 86)
@@ -408,5 +408,9 @@
 
 !Mars possible matter with dtphys in input and include!!!
+#ifdef MESOSCALE
+      CALL meso_inifis(1,llm,day0,daysec,dtphys,
+#else
       CALL inifis(1,llm,day0,daysec,dtphys,
+#endif
      .            lati,long,area,rad,g,r,cpp)
 c   Initialisation pour prendre en compte les vents en 1-D
@@ -591,6 +595,9 @@
 c       appel de la physique
 c       --------------------
-
+#ifdef MESOSCALE
+      CALL meso_physiq (1,llm,nqmx,
+#else
       CALL physiq (1,llm,nqmx,
+#endif
      ,     firstcall,lastcall,
      ,     day,time,dtphys,
Index: trunk/mesoscale/LMDZ.MARS.new/myGCM/DEFS_JB/callphys.def.orig
===================================================================
--- trunk/mesoscale/LMDZ.MARS.new/myGCM/DEFS_JB/callphys.def.orig	(revision 86)
+++ trunk/mesoscale/LMDZ.MARS.new/myGCM/DEFS_JB/callphys.def.orig	(revision 86)
@@ -0,0 +1,117 @@
+##General options
+##~~~~~~~~~~~~~~~
+#Run with or without tracer transport ?
+tracer=.true.
+
+#Diurnal cycle ?  if diurnal=False, diurnal averaged solar heating
+diurnal=.true.
+
+#Seasonal cycle ? if season=False, Ls stays constant, to value set in "start"
+season = .true.   
+
+#write some more output on the screen ?
+lwrite = .false.    
+
+#Save statistics in file "stats.nc" ?
+stats =.false.
+
+#Save EOF profiles in file "profiles" for Climate Database?
+calleofdump = .false.
+
+## Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+## ~~~~~~~~~~~~~
+#  =1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario,
+#  =4 Mars Year 24 from TES assimilation (same as =24 for now)
+#  =24 Mars Year 24 from TES assimilation (ie: MCD reference case)
+#  =25 Mars Year 25 from TES assimilation (ie: a year with a global dust storm)
+#  =26 Mars Year 26 from TES assimilation
+iaervar = 26 
+# Dust opacity at 700 Pa (when constant, i.e. for the iaervar=1 case)
+tauvis=0.2
+# Dust vertical distribution: 
+# (=0: old distrib. (Pollack90), =1: top set by "topdustref",
+#  =2: Viking scenario; =3 MGS scenario)
+iddist  = 3
+# Dust top altitude (km). (Matters only if iddist=1)
+topdustref = 55.
+
+## Physical Parameterizations :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call radiative transfer ?
+callrad   = .true.
+# call NLTE radiative schemes ?   matters only if callrad=T
+callnlte = .true.
+# call CO2 NIR absorption ?   matters only if callrad=T
+callnirco2 = .true.
+# call turbulent vertical diffusion ?
+calldifv  = .true.
+# call convective adjustment ?
+calladj  = .true.
+# call CO2 condensation ?
+callcond =.true.
+# call thermal conduction in the soil ?
+callsoil = .true.
+# call Lott's gravity wave/subgrid topography scheme ?
+calllott  = .true.
+
+## Radiative transfer options :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# the rad.transfer is computed every "iradia" physical timestep
+iradia    = 1
+# Output of the exchange coefficient mattrix ? for diagnostic only
+callg2d  = .false.
+# Rayleigh scattering : (should be .false. for now)
+rayleigh  = .false.
+
+## Tracer (dust water, ice and/or chemical species) options (used if tracer=T):
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# DUST: Transported dust ? (if >0, use 'dustbin' dust bins)
+dustbin    = 2
+# DUST: Radiatively active dust ? (matters if dustbin>0)
+active  = .true.
+# DUST: use mass and number mixing ratios to predict dust size ?
+doubleq   = .true.
+# DUST: use a small population of dust particules (submicron dust)?
+submicron = .false.
+# DUST: lifted by GCM surface winds ?
+lifting = .true.
+# DUST: lifted by dust devils ?
+callddevil = .false.
+# DUST: Scavenging by CO2 snowfall ?
+scavenging = .false.
+# DUST/WATERICE: Gravitationnal sedimentation ?
+sedimentation = .true.
+# WATERICE: Radiatively active transported atmospheric water ice ?
+activice   = .true.
+# WATER: Compute water cycle
+water = .true.
+# WATER: current permanent caps at both poles. True IS RECOMMENDED
+#        (with .true., North cap is a source of water and South pole
+#         is a cold trap)
+caps  = .true.
+# PHOTOCHEMISTRY: include chemical species 
+photochem  = .false.
+
+## Thermospheric options (relevant if tracer=T) :
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call thermosphere ?
+callthermos = .false.
+# WATER: included without cycle (only if water=.false.)
+thermoswater = .false.
+# call thermal conduction ? (only if callthermos=.true.)
+callconduct = .false.
+# call EUV heating ? (only if callthermos=.true.)
+calleuv=.false.
+# call molecular viscosity ? (only if callthermos=.true.)
+callmolvis = .false.
+# call molecular diffusion ? (only if callthermos=.true.)
+callmoldiff = .false.
+# call thermospheric photochemistry ? (only if callthermos=.true.)
+thermochem = .false.
+# date for solar flux calculation: (1985 < date < 2002)
+## (Solar min=1996.4 ave=1993.4 max=1990.6)
+solarcondate = 1993.4
+
+## albedoes of CO2 ice caps from TES ?
+TESicealbedo=.true.
+
Index: trunk/mesoscale/LMD_MM_MARS/SIMU/gnome_launch.meso
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SIMU/gnome_launch.meso	(revision 85)
+++ trunk/mesoscale/LMD_MM_MARS/SIMU/gnome_launch.meso	(revision 86)
@@ -13,4 +13,5 @@
 ## @ class =CP
 # @ environment=tmpdir=$(job_name).$(jobid).results
+# @ notification = never
 # @ queue
 grep "^# @" $0
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM	(revision 85)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM	(revision 86)
@@ -101,6 +101,7 @@
 state  real  HR_NIR    ikj  misc  1  -  rd   "HR_NIR"    "HEATING RATE nirco2"             "K/s"     #SAVEMARS3 zdtnirco2
 state  real  HR_NLTE   ikj  misc  1  -  rd   "HR_NLTE"   "HEATING RATE nlte"               "K/s"     #SAVEMARS3 zdtnlte
-state  real  ALBBARE   ij   misc  1  -  rhd  "ALBBARE"   "SOIL ALBEDO"                     ""        #SAVEMARS2 albedodat
+state  real  ALBBARE   ij   misc  1  -  rd   "ALBBARE"   "SOIL ALBEDO"                     ""        #SAVEMARS2 albedodat
 state  real  VMR_ICE   ikj  misc  1  -  rhd  "VMR_ICE"   "VOL. MIXING RATIO ICE"           "ppm"     #SAVEMARS3 vmr
+state  real  TAU_ICE   ij   misc  1  -  rhd  "TAU_ICE"   "CLOUD OD at 825 cm-1 TES"        ""        #SAVEMARS2 tauTES
 ####
 ####
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_g95
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_g95	(revision 85)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_g95	(revision 86)
@@ -181,6 +181,6 @@
 # Ehouarn 'g95' compiler
 # Aymeric : added option for NETCDF libraries
-   set optim90=" -O2 -fno-second-underscore -fstatic -fzero -Wall "
-   set optimtru90=" -O2 -fno-second-underscore -fstatic -fzero -Wall "
+   set optim90=" -DMESOSCALE -O2 -fno-second-underscore -fstatic -fzero -Wall "
+   set optimtru90=" -DMESOSCALE -O2 -fno-second-underscore -fstatic -fzero -Wall "
 #   set opt_link=" -Mfree -L/usr/local/pgi/linux86/lib -lpgf90 -lpgftnrtl -lpghpf -lpghpf2 -L$NCDFLIB -lnetcdf -Wl,-Bstatic "
 # Ehouarn: there is no /usr/local/pgi/linux86/lib, but it doesn't matter
@@ -378,7 +378,7 @@
 #           set optim90="-mcmodel=medium"
 #         Aymeric : thought I'd complete the g95 debugging 
-          set optim="$optim"" -g -fbounds-check "
-          set optim90="$optim90"" -g -fbounds-check "
-          set optimtru90="$optimtru90"" -g -fbounds-check "
+          set optim="$optim"" -g -fbounds-check -DMESOSCALE "
+          set optim90="$optim90"" -g -fbounds-check -DMESOSCALE "
+          set optimtru90="$optimtru90"" -g -fbounds-check -DMESOSCALE "
         else 
            echo "pas d option debug predefinie pour cette machine"
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_ifort
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_ifort	(revision 85)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_ifort	(revision 86)
@@ -181,7 +181,7 @@
 # Ehouarn 'ifort' compiler
    #NB: on gnome -O3 ==> NaNs ...
-   set optim=" -O2 -ip -mkl=sequential -align common "
-   set optim90=" -O2 -ip -mkl=sequential -align common "
-   set optimtru90=" -O2 -ip -mkl=sequential -align common "
+   set optim=" -DMESOSCALE -O2 -ip -mkl=sequential -align common "
+   set optim90=" -DMESOSCALE -O2 -ip -mkl=sequential -align common "
+   set optimtru90=" -DMESOSCALE -O2 -ip -mkl=sequential -align common "
 #   set opt_link=" -Mfree -lpgf90 -lpgftnrtl -lpghpf -lpghpf2 -L$NCDFLIB -lnetcdf -Bstatic "
    set mod_loc_dir=$LIBOGCM
@@ -374,7 +374,7 @@
         else if $LINUX then
            ## for ifort
-             set optim=" -g -no-ftz -traceback -ftrapuv -fp-stack-check "
-             set optim90=" -g -no-ftz -traceback -ftrapuv -fp-stack-check "
-             set optimtru90=" -g -no-ftz -traceback -ftrapuv -fp-stack-check "
+             set optim=" -DMESOSCALE -g -no-ftz -traceback -ftrapuv -fp-stack-check "
+             set optim90=" -DMESOSCALE -g -no-ftz -traceback -ftrapuv -fp-stack-check "
+             set optimtru90=" -DMESOSCALE -g -no-ftz -traceback -ftrapuv -fp-stack-check "
         else 
            echo "pas d option debug predefinie pour cette machine"
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpi
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpi	(revision 85)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpi	(revision 86)
@@ -189,10 +189,10 @@
 ##-Ktrap=fp
 #   set optim90=" -fast -I. -Mlarge_arrays -mcmodel=medium" 
-   set optim=" -O2 -I. -Munroll -Mcache_align" #-C #-fpic" #-Mlarge_arrays -mcmodel=medium
-   set optim90=" -O2 -I. -Munroll -Mcache_align" #-C #-fpic" #-Mlarge_arrays -mcmodel=medium
+   set optim=" -DMESOSCALE -O2 -I. -Munroll -Mcache_align" #-C #-fpic" #-Mlarge_arrays -mcmodel=medium
+   set optim90=" -DMESOSCALE -O2 -I. -Munroll -Mcache_align" #-C #-fpic" #-Mlarge_arrays -mcmodel=medium
 ## Aymeric: added handling of large arrays and large memory
 ## Aymeric: get rid of -Mfree for mpi
 #   set optimtru90=" -fast -c -Mlarge_arrays -mcmodel=medium" 
-  set optimtru90=" -O2 -I. -Munroll -Mcache_align" #-C #-fpic" #-Mlarge_arrays -mcmodel=medium
+  set optimtru90=" -DMESOSCALE -O2 -I. -Munroll -Mcache_align" #-C #-fpic" #-Mlarge_arrays -mcmodel=medium
 #   set opt_link=" -Mfree -L/usr/local/pgi/linux86/lib -lpgf90 -lpgftnrtl -lpghpf -lpghpf2 -L$NCDFLIB -lnetcdf -Wl,-Bstatic "
 # Ehouarn: there is no /usr/local/pgi/linux86/lib, but it doesn't matter
@@ -388,5 +388,5 @@
            set optim90="$optim90"" -G1 "
         else if $LINUX then
-           set optim90="-g -Mbounds "
+           set optim90="-DMESOSCALE -g -Mbounds "
         else 
            echo "pas d option debug predefinie pour cette machine"
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpifort
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpifort	(revision 85)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpifort	(revision 86)
@@ -181,7 +181,7 @@
 # Ehouarn 'ifort' compiler
    #NB: on gnome -O3 ==> NaNs ...
-   set optim=" -O2 -ip -mkl=sequential -align common "
-   set optim90=" -O2 -ip -mkl=sequential -align common "
-   set optimtru90=" -O2 -ip -mkl=sequential -align common "
+   set optim=" -DMESOSCALE -O2 -ip -mkl=sequential -align common "
+   set optim90=" -DMESOSCALE -O2 -ip -mkl=sequential -align common "
+   set optimtru90=" -DMESOSCALE -O2 -ip -mkl=sequential -align common "
 #   set opt_link=" -Mfree -lpgf90 -lpgftnrtl -lpghpf -lpghpf2 -L$NCDFLIB -lnetcdf -Bstatic "
    set mod_loc_dir=$LIBOGCM
@@ -374,7 +374,7 @@
         else if $LINUX then
            ## for ifort
-             set optim=" -g -no-ftz -traceback -ftrapuv -fp-stack-check "
-             set optim90=" -g -no-ftz -traceback -ftrapuv -fp-stack-check "
-             set optimtru90=" -g -no-ftz -traceback -ftrapuv -fp-stack-check "
+             set optim=" -DMESOSCALE -g -no-ftz -traceback -ftrapuv -fp-stack-check "
+             set optim90=" -DMESOSCALE -g -no-ftz -traceback -ftrapuv -fp-stack-check "
+             set optimtru90=" -DMESOSCALE -g -no-ftz -traceback -ftrapuv -fp-stack-check "
         else 
            echo "pas d option debug predefinie pour cette machine"
Index: trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_pgf
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_pgf	(revision 85)
+++ trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_pgf	(revision 86)
@@ -180,6 +180,6 @@
 else if $LINUX then
 ##-Ktrap=fp
-   set optim90=" -fast -I." 
-   set optimtru90=" -fast -c -Mfree" 
+   set optim90=" -DMESOSCALE -fast -I." 
+   set optimtru90=" -DMESOSCALE -fast -c -Mfree" 
 #   set opt_link=" -Mfree -L/usr/local/pgi/linux86/lib -lpgf90 -lpgftnrtl -lpghpf -lpghpf2 -L$NCDFLIB -lnetcdf -Wl,-Bstatic "
 # Ehouarn: there is no /usr/local/pgi/linux86/lib, but it doesn't matter
@@ -379,5 +379,5 @@
            set optim90="$optim90"" -G1 "
         else if $LINUX then
-           set optim90="-g -Mbounds "
+           set optim90="-g -Mbounds -DMESOSCALE "
         else 
            echo "pas d option debug predefinie pour cette machine"
Index: trunk/mesoscale/LMD_MM_MARS/makemeso
===================================================================
--- trunk/mesoscale/LMD_MM_MARS/makemeso	(revision 85)
+++ trunk/mesoscale/LMD_MM_MARS/makemeso	(revision 86)
@@ -439,6 +439,6 @@
        touch callkeys.h.meso  
          ln -sf meso_callkeys.h callkeys.h
-       touch dustlift.F.meso 
-         ln -sf meso_dustlift.F dustlift.F  ## attention avec "ancienne nouvelle physique"
+       #touch dustlift.F.meso 
+       #  ln -sf meso_dustlift.F dustlift.F  ## attention avec "ancienne nouvelle physique"
        touch readtesassim.F90.meso
          ln -sf meso_readtesassim.F90 readtesassim.F90  ## attention idem
@@ -467,7 +467,7 @@
   then
     echo 'DEBUG DEBUG DEBUG DEBUG'
-    nohup makegcm_${compilo} -debug -t ${tra} -p mars -d ${physz} meso_testphys1d | tee libo/log_compile_phys | grep 'warnings' >> libo/log_compile_phys 
+    nohup makegcm_${compilo} -debug -t ${tra} -p mars -d ${physz} testphys1d | tee libo/log_compile_phys | grep 'warnings' >> libo/log_compile_phys 
   else
-    nohup makegcm_${compilo} -t ${tra} -p mars -d ${physz} meso_testphys1d | tee libo/log_compile_phys | grep 'warnings' >> libo/log_compile_phys
+    nohup makegcm_${compilo} -t ${tra} -p mars -d ${physz} testphys1d | tee libo/log_compile_phys | grep 'warnings' >> libo/log_compile_phys
   fi
   echo ... done
Index: trunk/mesoscale/NOTES.txt
===================================================================
--- trunk/mesoscale/NOTES.txt	(revision 85)
+++ trunk/mesoscale/NOTES.txt	(revision 86)
@@ -1,2 +1,5 @@
+
+A FAIRE:::: mettre des flags precompilo dans les meso_
+les reporter dans makegcm
 
 changer le renormalisation dans aeropacity ????
Index: trunk/mesoscale/PLOT/MINIMAL/map_latlon.pro
===================================================================
--- trunk/mesoscale/PLOT/MINIMAL/map_latlon.pro	(revision 85)
+++ trunk/mesoscale/PLOT/MINIMAL/map_latlon.pro	(revision 86)
@@ -389,5 +389,5 @@
 	c_labels=findgen(n_elements(lev))*0., $
 ;        c_charsize=!P.charsize/2, $
-	color=0;50;75;100tropclair;150;0;, $
+	color=75;0;50;75;100tropclair;150;0;, $
 ;	/noerase, $
 ;	xtitle=xtitleset, $
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/movie.idl
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/movie.idl	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/movie.idl	(revision 86)
@@ -0,0 +1,85 @@
+;;---------------------------------------
+;  MOVIE
+;
+;  **** DO NOT MODIFY ****
+;
+;  set parameters in plot.def then type:
+;
+;  @movie.idl
+;
+;;---------------------------------------
+
+field1='' & field2='' & winds=''
+backup_data='no' & already_data='no'  	
+datafolder='./' & plotfolder='./'
+nomovie='yes'
+topo=0
+extract='no'
+
+SPAWN, '\rm user.idl'
+SPAWN, 'grep -v grep plot.def > user.idl'
+
+;; to avoid dumb and useless warnings from IDL
+space=0. & altitude=0. & lon=0. & lat=0.
+what_I_plot=0. & title_user='' & missing_value=0.
+
+@user.idl
+if (nam5 eq " interp_method = 0 ") then coord='height'
+if ((nam5 eq " interp_method = 1 ") and (tabnam(0) lt tabnam(1))) then coord='height'
+
+.compile call_arwpost
+if (extract eq 'yes') then call_arwpost, nam1, nam2, nam3, nam4, nam5, tabnam
+
+denom=plot+'_'+coord+'_'+field1
+if ((field2 ne '') and (topo eq 0)) then denom=denom+'_'+field2
+if (topo eq 1) then denom=denom+'_HGT'
+if (winds(0) ne '') then denom=denom+'_'+winds(0)+winds(1)
+
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.ctl input.ctl'
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.dat input.dat'
+
+SPAWN, '\rm param_plot.idl'
+SPAWN, 'grep grepALL plot.def > param_plot.idl'
+if ((plot eq 'meridional') or (plot eq 'zonal')) then SPAWN, 'grep grepSEC plot.def >> param_plot.idl'
+if (plot eq 'map') then SPAWN, 'grep grepMAP plot.def >> param_plot.idl'
+.compile map_latlon
+.compile section
+
+!p.charthick = 2.0
+!p.thick = 3.0
+!x.thick = 2.0
+!y.thick = 2.0
+
+for i=1,num do $
+        out_wrf, $
+        plot=plot, $
+        field1=field1, $
+        field2=field2, $
+        when=i, $
+        winds=winds, $
+        level=level,$
+        save_data=yeah,$        ;; data is loaded only once
+        topo=topo,$
+;        range=range,$
+        colors=colors, $
+        nlon=nlon, $
+        nlat=nlat, $
+        save_ps=denom+string(1000+i,'(I0)')
+
+print, 'generating movie ...'
+SPAWN, 'rm -rf temp'
+SPAWN, 'mkdir temp'
+SPAWN, 'mv -f '+denom+'????.ps temp'
+SPAWN, 'cp bigconvert200 temp'
+SPAWN, 'cd temp ; bigconvert200 *.ps ; convert -delay 60 *.png movie.gif ; cd ..'
+if (nomovie ne 'yes') then SPAWN, 'animate temp/movie.gif &' 
+print, 'movie done'
+
+SPAWN, 'mv -f temp/'+denom+'????.ps '+plotfolder+'/'
+SPAWN, 'mv -f temp/movie.gif '+plotfolder+'/'+denom+'.gif'
+SPAWN, 'mv -f temp/'+denom+'????.png '+plotfolder+'/'
+if (backup_data eq 'yes') then SPAWN, 'cp -f namelist.ARWpost '+datafolder+'/'+coord+'.namelist.ARWpost'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.ctl '+datafolder+'/'+coord+'.input.ctl'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.dat '+datafolder+'/'+coord+'.input.dat'
+
+SPAWN, 'cp -f plot.def '+plotfolder+'/'+denom+'.def'
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/out_diagfi.pro
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/out_diagfi.pro	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/out_diagfi.pro	(revision 86)
@@ -0,0 +1,471 @@
+pro out_diagfi,	$
+		plot=plot_user, $		; 'map', 'profile' (default is 'map')
+	;----------------
+	; input fields
+	;----------------
+		field1=field_user, $		; field to contour: 'u', 'v', 'temp', ... etc ...
+		field2=fieldc_user, $		; field to overplot (contour lines, optional)
+	;------------------------------------------------------
+	; 'map' options (no effect when plot is 'profile')
+	;------------------------------------------------------
+		topo=topo, $			; overplot topography contours (!! must load phisinit !! override field2)
+		winds=winds, $			; option to overplot winds (e.g. =['U','V'])
+		level=level, $			; vertical level subscript (starts at 1 ...)(if -1: all levels)
+	;------------------------------------------------------
+	; 'profile' options (no effect when plot is 'map')
+	;------------------------------------------------------
+		inprofile=inprofile, $		; a profile to compare with ... (override field2)
+		incolumn=incolumn, $		; altitudes ... in case /= column
+		;; NB: inprofile and incolumn may be used to output profiles
+		;; ... reset with inprofile=0, incolumn=0
+	;----------------
+	; choose time
+	;----------------
+		when=when, $			; time subscript (starts at 1 ...)(if -1: first time step)
+	;----------------
+	; plot options
+	;----------------
+		range=range_user, $		; range=[min,max]
+		colors=colors, $		; 16, 32, 64, ... etc
+	;-------------------------
+	; no useless data loading
+	;-------------------------
+		save_data=save_data		; to avoid loading data if already done once
+
+
+;---------------------------------------------------------------------------
+;
+; A general routine to plot diagfi.nc GCM outputs
+;
+; USE: see out_wrf
+;
+;---------------------------------------------------------------------------
+; A. Spiga, August 2007, adding plot keyword September 2007
+;---------------------------------------------------------------------------
+
+
+
+;-------------------------------------
+; default settings and user settings
+;-------------------------------------
+charend=''
+diff_scale=0  ;1
+if (n_elements(plot_user) eq 0) then plot_user='map'
+if (n_elements(field_user) eq 0) then field_user='all'
+if (n_elements(level) eq 0) then level=-1
+if (n_elements(when) eq 0) then when=-1
+if (n_elements(colors) eq 0) then colors=32
+if (n_elements(minfield_init)*n_elements(maxfield_init) eq 0) then begin
+	minfield_init = 0.
+	maxfield_init = 0.
+endif
+
+if (n_elements(range_user) eq 2) then begin
+	minfield_init=float(range_user(0))
+	maxfield_init=float(range_user(1))
+	print, 'setting user limits !', minfield_init, maxfield_init 
+	diff_scale=1
+endif
+
+
+
+;******************************************************************************
+;
+; BEGIN: different from out_wrf
+;
+;******************************************************************************
+
+;TEMP TEMP
+vert_coord='sigma'
+;vert_coord='pressure'
+;vert_coord='height'
+;TEMP TEMP
+
+;----------------------
+; set parameters
+;----------------------
+file='input.diagfi'
+if (n_elements(topo) ne 0) then fieldc_user='phisinit'	;case with topography
+
+
+
+;----------------------------------------------------
+; Open diagfi file GCM and get coordinates & field
+;----------------------------------------------------
+cdfid = ncdf_open(file)
+
+print, 'read coord ...'
+varid=ncdf_varid(cdfid,'latitude')
+        ncdf_varget, cdfid, varid, lat_gcm
+varid=ncdf_varid(cdfid,'longitude')
+        ncdf_varget, cdfid, varid, lon_gcm
+varid=ncdf_varid(cdfid,'Time')
+        ncdf_varget, cdfid, varid, hour
+lon=lon_gcm
+west_east=n_elements(lon)
+lat=lat_gcm
+south_north=n_elements(lat)
+print, 'done !'
+;;print, 'time ...', hour*24 MOD 24, ' hours and day number ', ceil(hour)
+
+; field
+if (n_elements(save_data) eq 0) then begin
+	print, 'read data ...'
+	varid=ncdf_varid(cdfid,field_user)
+		ncdf_varget, cdfid, varid, var
+		bottom_top=n_elements(var(0,0,*,0))
+		time=n_elements(var(0,0,0,*))
+;save_array=fltarr(3,west_east,south_north,bottom_top,time)
+;save_array(0,*,*,*,*)=var(*,*,*,*)
+;;**NB: ralentit beaucoup
+	print, 'done !'
+endif else begin	
+	var=reform(save_data(0,*,*,*,*))
+		bottom_top=n_elements(var(0,0,*,0))
+		time=n_elements(var(0,0,0,*))
+endelse
+
+print, 'vertical levels in file ... ', bottom_top
+print, 'time steps in file ...', time
+
+
+
+	;-----------------------------
+	; min/max on the whole field
+	;-----------------------------
+	w=where(abs(var) lt 1e30)
+	if (w(0) ne -1) then begin
+		minfield=min(var[w])
+		maxfield=max(var[w])
+			;;; manage unvalid points
+			;;;w=where(var eq 1e37)
+			;;;if (w(0) ne -1) then var[w]=!Values.F_NAN
+	endif else begin
+		var=0
+		print, 'no valid values in this field', charvar 
+	endelse
+
+
+
+;---------------------
+; Get altitude
+;---------------------
+if (vert_coord ne 'pressure') then begin
+	print, 'read altitude ...'
+	varid=ncdf_varid(cdfid,'aps')
+        	ncdf_varget, cdfid, varid, aps
+	varid=ncdf_varid(cdfid,'bps')
+        	ncdf_varget, cdfid, varid, bps
+	varid=ncdf_varid(cdfid,'ps')
+        	ncdf_varget, cdfid, varid, ps
+	varid=ncdf_varid(cdfid,'altitude')
+        	ncdf_varget, cdfid, varid, z
+	print, 'done !'
+endif else begin
+	print, 'read altitude ...'
+	varid=ncdf_varid(cdfid,'altitude')
+        	ncdf_varget, cdfid, varid, z
+	z=z/100.   ;mbar
+	print, 'done !'
+endelse
+
+
+;-------------
+; Get winds
+;-------------
+if (n_elements(winds) ne 0) then begin
+if (n_elements(save_data) eq 0) then begin
+	print, 'read winds ...'
+	varid=ncdf_varid(cdfid,winds(0))
+		ncdf_varget, cdfid, varid, u
+	varid=ncdf_varid(cdfid,winds(1))
+		ncdf_varget, cdfid, varid, v
+	;save_array(1,*,*,*,*)=u(*,*,*,*)
+	;save_array(2,*,*,*,*)=v(*,*,*,*)
+	print, 'done !'
+endif else begin	
+	u=reform(save_data(1,*,*,*,*))
+	v=reform(save_data(2,*,*,*,*))
+endelse
+endif
+
+;-------------------------------
+; Get another field to contour
+;-------------------------------
+if (n_elements(fieldc_user) ne 0) then begin
+	print, 'get alternate field ...'
+	varid=ncdf_varid(cdfid,fieldc_user)
+		ncdf_varget, cdfid, varid, varcontour
+	print, 'done !'
+endif
+
+
+;;-------------------------------
+;; Save data for further use ...
+;;-------------------------------
+;if (n_elements(save_data) eq 0) then begin
+;	save_data=save_array
+;	save_array=0. ;free 'save_array' array
+;endif
+
+
+;---------------------------------
+; Some stuff for continuity
+;---------------------------------
+;;if (vert_coord ne 'pressure') then surfpres=reform(ps(*,*,nt))
+;;A CORRIGER - A METTRE DANS LA BOUCLE
+if (vert_coord ne 'pressure') then surfpres=reform(ps(*,*,0))
+if (vert_coord eq 'sigma') then begin
+	z=aps/surfpres(0,0)+bps
+	z=z*100.
+	vert_coord='model_level'
+endif
+charend='_diagfi'
+
+case field_user of 
+'temp': charvar='tk'
+'u': charvar='U'
+'v': charvar='V'
+'ps': charvar='PSFC'
+else:
+endcase
+
+;;
+;;
+
+
+;******************************************************************************
+;
+; END: different from out_wrf
+;
+;******************************************************************************
+
+
+
+;----------------------
+; initialize graphics
+;----------------------
+
+!p.charthick = 2.0
+!p.thick = 3.0
+!x.thick = 2.0
+!y.thick = 2.0
+
+if (when ne -1) then charend=charend+string(when, '(I0)')
+
+set_plot, 'ps'
+device, filename=plot_user+'_'+vert_coord+'_'+charvar+charend+'.ps', /color
+
+case charvar of
+'tk': begin
+	title='Temperature (K)'
+	pal=3
+end
+'U': begin
+	title='Zonal wind (!Nm!N.s!U-1!N)'
+	pal=33
+end
+'V': begin
+	title='Zonal wind (!Nm!N.s!U-1!N)'
+	pal=33
+end
+'pressure': begin
+	title='Pressure (hPa)'
+	pal=16
+end
+'height': begin
+	title='Height (km)'
+	pal=16
+end
+'T': title='Potential temperature departure from t0 (K)'
+'W': title='Vertical wind (!Ncm!N.s!U-1!N)'
+else: title=''
+endcase
+title_save=title
+
+
+
+;---------------------------------
+; choose positioning ...
+;---------------------------------
+
+case plot_user of
+'map': begin
+;
+; 1. map: looping on vertical levels
+;
+	if (level eq -1) then begin
+		theloopmin=0
+		theloopmax=bottom_top-1
+	endif else begin
+		theloopmin=level-1 ;!! IDL arrays starts at 0
+		theloopmax=level-1
+	endelse	
+	if (when eq -1) then nt = 0 else nt=when-1
+end
+'profile': begin
+;
+; 2. profile: looping on time
+;
+	;;discrete=4
+	title_axis=strarr(2)
+	title_axis(0)=title_save
+	if (when eq -1) then begin
+		theloopmin=0
+		theloopmax=time-1
+	endif else begin
+		theloopmin=when-1 ;!! IDL arrays starts at 0
+		theloopmax=when-1
+	endelse	
+
+	column=z
+	case vert_coord of
+		'height': title_axis(1)='Altitude (km)'
+		'pressure': begin
+			H=10.
+			p0=6.1
+			column=H*alog(p0/z)
+			title_axis(1)='Log-Pressure altitude (km)'
+		end
+		'model_level': title_axis(1)='Model level'
+	endcase
+
+;****TEMP
+nx=round(west_east/2)
+ny=round(south_north/2)
+;****TEMP
+
+	displaylat=float(round(lat(ny)*100))/100.
+	displaylon=float(round(lon(nx)*100))/100.
+	description='WRF Profile at latitude '+string(displaylat,format='(1F6.2)')+$
+	'!Uo'+'!N and longitude '+string(displaylon,format='(1F7.2)')+'!Uo!N '
+	title=description
+
+end
+else: stop
+endcase
+
+
+;-------------
+; plot loop
+;-------------	
+for theloop=theloopmin,theloopmax do begin
+
+case plot_user of
+'map': begin
+	;
+	; case 1: map
+	;	
+	case vert_coord of
+	'height': charlev=string(round(z(theloop)*1000.),'(I0)')+' m'
+	'pressure': charlev=string(round(z(theloop)*100.),'(I0)')+' Pa'
+	'model_level': charlev=string(round(z(theloop)),'(I0)')
+	endcase
+	title=title_save+' at level '+charlev
+	print, title
+
+	;-------------
+	; draw map
+	;-------------
+
+	what_I_plot=var[*,*,theloop,nt]
+
+	if (diff_scale eq 0) then begin
+		minfield_init=minfield	;-same scale for all levels 
+		maxfield_init=maxfield	;-same scale for all levels
+	endif
+	if (n_elements(fieldc_user) eq 0) then begin
+		overcontour=0.
+	endif else begin
+		if (n_elements(topo) ne 0) then begin
+			overcontour=varcontour[*,*,0,0] ; other levels than 0 are blank !
+		endif else begin
+			overcontour=varcontour[*,*,theloop,nt]
+		endelse
+	endelse
+	if (n_elements(winds) eq 0) then begin
+		overvector_x=0.
+		overvector_y=0.
+	endif else begin
+		overvector_x=u[*,*,theloop,nt]
+		overvector_y=v[*,*,theloop,nt]
+	endelse
+
+
+	map_latlon, $
+		what_I_plot, $
+		lon, $
+		lat, $
+		ct=pal, $
+		minfield=minfield_init, $
+		maxfield=maxfield_init, $
+		overcontour=overcontour, $
+		overvector_x=overvector_x, $
+		overvector_y=overvector_y, $
+		colors=colors, $
+		title=title
+end
+'profile': begin
+	;
+	; case 2: profile
+	;
+	print, title
+
+	what_I_plot=var[nx,ny,*,theloop]
+	if (diff_scale eq 0) then begin
+		minfield_init=minfield	;-same scale for all levels 
+		maxfield_init=maxfield	;-same scale for all levels
+	endif
+	if (n_elements(fieldc_user) eq 0) then begin
+		overcontour=0.
+		overcontourz=0.
+	endif else begin
+		overcontour=varcontour[nx,ny,*,theloop]
+		overcontourz=0.
+	endelse
+	if (n_elements(inprofile) gt 1) then begin
+		overcontour=inprofile
+		overcontourz=0.
+		if (n_elements(incolumn) gt 1) then overcontourz=incolumn
+	endif
+
+	profile, $
+		what_I_plot, $				; 1D vertical profile
+		column, $				; altitudes	
+	;	alt=alt, $				; altitude range [altmin, altmax]
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		inprofile=overcontour, $		; another vertical profile to overplot
+		incolumn=overcontourz, $		; altitudes of the other vertical profile (in case /= column)
+		title_plot=title, $			; title of the plot ('Profile' is default)
+		title_axis=title_axis, $		; title of the [x,y] axis (['Field','Altitude'] is default)
+	;	mention=mention, $			; add text precision within the plot window (default is nothing or '')
+		discrete=discrete			; show the profile points (= type of points in !psym)
+end
+else: 
+endcase
+
+
+
+endfor
+device, /close
+
+	; save the profiles for another plot (last of the time loop)
+	if ((n_elements(inprofile) le 1) and (plot_user eq 'profile')) then begin
+		inprofile=what_I_plot
+		incolumn=column
+	endif
+
+
+;******************************************************************************
+;
+; BEGIN: different from out_wrf
+;
+;******************************************************************************
+;;;skip:
+;;;endfor
+
+end
+
+
+
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/out_diagfi_new.pro
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/out_diagfi_new.pro	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/out_diagfi_new.pro	(revision 86)
@@ -0,0 +1,768 @@
+pro out_diagfi, $
+		plot=plot_user, $		; 'map', 'profile', 'hovmuller', 'meridional', 'zonal' (default is 'map')
+	;----------------
+	; input fields
+	;----------------
+		field1=field_user, $		; field to contour: 'U', 'V', 'tk', ... etc ...
+		field2=fieldc_user, $		; field to overplot (contour lines, optional)
+	;------------------------------------------------------
+	; 'map' options (no effect when plot is 'profile')
+	;------------------------------------------------------
+		topo=topo, $			; overplot topography contours (!! must load HGT !! override field2)
+		winds=winds, $			; option to overplot winds (e.g. =['U','V'])
+		level=level, $			; vertical level subscript (starts at 1 ...)(if -1: all levels)
+	;------------------------------------------------------
+	; 'profile' options (no effect when plot is 'map')
+	;------------------------------------------------------
+		inprofile=inprofile, $		; a profile to compare with ... (override field2)
+		incolumn=incolumn, $		; altitudes ... in case /= column
+		;; NB: inprofile and incolumn may be used to output profiles
+		;; ... reset with inprofile=0, incolumn=0
+	;----------------
+	; choose time
+	;----------------
+		when=when, $			; time subscript (starts at 1 ...)(if -1: first time step)
+	;----------------
+	; choose latitude/longitude
+	;----------------
+		nlat=nlat, $			; latitude subscript (starts at 1 ...)(if -1: all latitudes)(default is middle)
+		nlon=nlon, $			; longitude subscript (starts at 1 ...)(if -1: all longitudes)(default is middle)
+	;----------------
+	; plot options
+	;----------------
+		range=range_user, $		; range=[min,max]
+		colors=colors, $		; 16, 32, 64, ... etc
+	;-------------------------
+	; no useless data loading
+	;-------------------------
+		save_data=save_data		; to avoid loading data if already done once
+
+
+
+;---------------------------------------------------------------------------
+;
+; A general routine to plot ARWpost outputs
+;
+; USE:
+; - out_wrf	(will plot everything !)
+; - out_wrf, field1='tk', range=[160,300]
+;
+; - out_wrf, field1='pressure', field2='tk'
+; - out_wrf, field1=TSK, field2=HGT, winds=['U','V'] 
+;
+; - out_wrf, save_data=field
+;   ... then next calls including 'save_data=field' won't include the data loading 
+;
+; - out_wrf, field1='tk', winds=['umet','vmet'], topo=1, save_data=yeah, level=1, when=2
+;
+; - out_wrf, field1='tk', plot='profile'
+;
+; - out_wrf, plot='profile', field1='tk', when=1, inprofile=inprofile, incolumn=incolumn
+;   out_diagfi, plot='profile', field1='temp', when=112, inprofile=inprofile, incolumn=incolumn 	
+;
+; - out_wrf, plot='meridional',field1='U'
+;
+; - out_wrf, plot='hovmuller',field1='TSK'
+;
+; - out_wrf, plot='zonal',field1='V', nlat=-1
+; 
+;---------------------------------------------------------------------------
+; A. Spiga, August 2007, adding plot keyword September 2007
+;---------------------------------------------------------------------------
+
+
+    
+
+
+;-------------------------------------
+; default settings and user settings
+;-------------------------------------
+diff_scale=0  ;;1
+charend=''
+format=''
+if (n_elements(plot_user) eq 0) then plot_user='map'
+if (n_elements(field_user) eq 0) then field_user='all'
+if (n_elements(level) eq 0) then level=-1
+if (n_elements(when) eq 0) then when=-1
+;;if (n_elements(nlon) eq 0) then nlon=-1
+;;if (n_elements(nlat) eq 0) then nlat=-1
+if (n_elements(colors) eq 0) then colors=32
+if (n_elements(minfield_init)*n_elements(maxfield_init) eq 0) then begin
+	minfield_init = 0.
+	maxfield_init = 0.
+endif
+
+if (n_elements(range_user) eq 2) then begin
+	minfield_init=float(range_user(0))
+	maxfield_init=float(range_user(1))
+	print, 'setting user limits !', minfield_init, maxfield_init 
+	diff_scale=1
+endif
+
+if (n_elements(topo) ne 0) then fieldc_user='HGT'	;case with topography
+
+
+
+;******************************************************************************
+;
+; BEGIN: different from out_wrf
+;
+;******************************************************************************
+
+;TEMP TEMP
+vert_coord='sigma'
+;vert_coord='pressure'
+;vert_coord='height'
+;TEMP TEMP
+
+;----------------------
+; set parameters
+;----------------------
+file='input.diagfi'
+if (n_elements(topo) ne 0) then fieldc_user='phisinit'	;case with topography
+
+
+
+;----------------------------------------------------
+; Open diagfi file GCM and get coordinates & field
+;----------------------------------------------------
+cdfid = ncdf_open(file)
+
+print, 'read coord ...'
+varid=ncdf_varid(cdfid,'latitude')
+        ncdf_varget, cdfid, varid, lat_gcm
+varid=ncdf_varid(cdfid,'longitude')
+        ncdf_varget, cdfid, varid, lon_gcm
+varid=ncdf_varid(cdfid,'Time')
+        ncdf_varget, cdfid, varid, hour
+lon=lon_gcm
+west_east=n_elements(lon)
+lat=lat_gcm
+south_north=n_elements(lat)
+print, 'done !'
+;;print, 'time ...', hour*24 MOD 24, ' hours and day number ', ceil(hour)
+
+; field
+if (n_elements(save_data) eq 0) then begin
+	print, 'read data ...'
+	varid=ncdf_varid(cdfid,field_user)
+		ncdf_varget, cdfid, varid, var
+		bottom_top=n_elements(var(0,0,*,0))
+		time=n_elements(var(0,0,0,*))
+;save_array=fltarr(3,west_east,south_north,bottom_top,time)
+;save_array(0,*,*,*,*)=var(*,*,*,*)
+;;**NB: ralentit beaucoup
+	print, 'done !'
+endif else begin	
+	var=reform(save_data(0,*,*,*,*))
+		bottom_top=n_elements(var(0,0,*,0))
+		time=n_elements(var(0,0,0,*))
+endelse
+
+print, 'vertical levels in file ... ', bottom_top
+print, 'time steps in file ...', time
+
+
+
+	;-----------------------------
+	; min/max on the whole field
+	;-----------------------------
+	w=where(abs(var) lt 1e30)
+	if (w(0) ne -1) then begin
+		minfield=min(var[w])
+		maxfield=max(var[w])
+			;;; manage unvalid points
+			;;;w=where(var eq 1e37)
+			;;;if (w(0) ne -1) then var[w]=!Values.F_NAN
+	endif else begin
+		var=0
+		print, 'no valid values in this field', charvar 
+	endelse
+
+
+
+;---------------------
+; Get altitude
+;---------------------
+if (vert_coord ne 'pressure') then begin
+	print, 'read altitude ...'
+	varid=ncdf_varid(cdfid,'aps')
+        	ncdf_varget, cdfid, varid, aps
+	varid=ncdf_varid(cdfid,'bps')
+        	ncdf_varget, cdfid, varid, bps
+	varid=ncdf_varid(cdfid,'ps')
+        	ncdf_varget, cdfid, varid, ps
+	varid=ncdf_varid(cdfid,'altitude')
+        	ncdf_varget, cdfid, varid, z
+	print, 'done !'
+endif else begin
+	print, 'read altitude ...'
+	varid=ncdf_varid(cdfid,'altitude')
+        	ncdf_varget, cdfid, varid, z
+	z=z/100.   ;mbar
+	print, 'done !'
+endelse
+
+
+;-------------
+; Get winds
+;-------------
+if (n_elements(winds) ne 0) then begin
+if (n_elements(save_data) eq 0) then begin
+	print, 'read winds ...'
+	varid=ncdf_varid(cdfid,winds(0))
+		ncdf_varget, cdfid, varid, u
+	varid=ncdf_varid(cdfid,winds(1))
+		ncdf_varget, cdfid, varid, v
+	;save_array(1,*,*,*,*)=u(*,*,*,*)
+	;save_array(2,*,*,*,*)=v(*,*,*,*)
+	print, 'done !'
+endif else begin	
+	u=reform(save_data(1,*,*,*,*))
+	v=reform(save_data(2,*,*,*,*))
+endelse
+endif
+
+;-------------------------------
+; Get another field to contour
+;-------------------------------
+if (n_elements(fieldc_user) ne 0) then begin
+	print, 'get alternate field ...'
+	varid=ncdf_varid(cdfid,fieldc_user)
+		ncdf_varget, cdfid, varid, varcontour
+	print, 'done !'
+endif
+
+
+;;-------------------------------
+;; Save data for further use ...
+;;-------------------------------
+;if (n_elements(save_data) eq 0) then begin
+;	save_data=save_array
+;	save_array=0. ;free 'save_array' array
+;endif
+
+
+;---------------------------------
+; Some stuff for continuity
+;---------------------------------
+;;if (vert_coord ne 'pressure') then surfpres=reform(ps(*,*,nt))
+;;A CORRIGER - A METTRE DANS LA BOUCLE
+if (vert_coord ne 'pressure') then surfpres=reform(ps(*,*,0))
+if (vert_coord eq 'sigma') then begin
+	z=aps/surfpres(0,0)+bps
+	z=z*100.
+	vert_coord='model_level'
+endif
+charend='_diagfi'
+
+case field_user of 
+'temp': charvar='tk'
+'u': charvar='U'
+'v': charvar='V'
+'ps': charvar='PSFC'
+else:
+endcase
+
+;;
+;;
+
+
+;******************************************************************************
+;
+; END: different from out_wrf
+;
+;******************************************************************************
+
+
+
+;----------------------
+; initialize graphics
+;----------------------
+
+!p.charthick = 2.0
+!p.thick = 3.0
+!x.thick = 2.0
+!y.thick = 2.0
+
+if (when ne -1) then charend=charend+string(when, '(I0)')
+
+set_plot, 'ps'
+device, filename=plot_user+'_'+vert_coord+'_'+charvar+charend+'.ps', /color
+
+case charvar of
+'HGT': begin
+        title='Topography (km)'
+	pal=16
+	level=1
+	format='(F7.0)'
+end
+'tk': begin
+	title='Temperature (K)'
+	pal=3
+	format='(F4.0)'
+end
+'TSK': begin
+	title='Surface Temperature (K)'
+	pal=3
+	level=1
+end
+'U': begin
+	title='Zonal wind (!Nm!N.s!U-1!N)'
+	pal=33
+end
+'V': begin
+	title='Meridional wind (!Nm!N.s!U-1!N)'
+	pal=33
+end
+'pressure': begin
+	title='Pressure (hPa)'
+	pal=16
+end
+'height': begin
+	title='Height (km)'
+	pal=16
+end
+'T': title='Potential temperature departure from t0 (K)'
+'W': begin
+	title='Vertical wind (!Nm!N.s!U-1!N)'
+	pal=16
+	format='(F6.2)'
+end
+else: begin
+	title=''
+	pal=33
+end
+endcase
+title_save=title
+
+
+
+;---------------------------------
+; choose positioning ...
+;---------------------------------
+
+case plot_user of
+'map': begin
+;
+; 1. map: looping on vertical levels
+;
+	if (level eq -1) then begin
+		theloopmin=0
+		theloopmax=bottom_top-1
+	endif else begin
+		theloopmin=level-1 ;!! IDL arrays starts at 0
+		theloopmax=level-1
+	endelse	
+	if (when eq -1) then nt = 0 else nt=when-1
+end
+'profile': begin
+;
+; 2. profile: looping on time
+;
+	;;discrete=4
+	title_axis=strarr(2)
+	title_axis(0)=title_save
+	if (when eq -1) then begin
+		theloopmin=0
+		theloopmax=time-1
+	endif else begin
+		theloopmin=when-1 ;!! IDL arrays starts at 0
+		theloopmax=when-1
+	endelse	
+
+	column=z
+	case vert_coord of
+		'height': title_axis(1)='Altitude (km)'
+		'pressure': begin
+			H=10.
+			p0=6.1
+			column=H*alog(p0/z)
+			title_axis(1)='Log-Pressure altitude (km)'
+		end
+		'model_level': title_axis(1)='Model level'
+	endcase
+
+	displaylat=float(round(lat(nlat)*100))/100.
+	displaylon=float(round(lon(nlon)*100))/100.
+	description='WRF Profile at latitude '+string(displaylat,format='(1F6.2)')+$
+	'!Uo'+'!N and longitude '+string(displaylon,format='(1F7.2)')+'!Uo!N '
+	title=description
+end
+'hovmuller': begin
+;
+; 3. hovmuller: looping on vertical levels
+;
+	if (level eq -1) then begin
+		theloopmin=0
+		theloopmax=bottom_top-1
+	endif else begin
+		theloopmin=level-1 ;!! IDL arrays starts at 0
+		theloopmax=level-1
+	endelse	
+end
+'zonal': begin
+;
+; 4. zonal: looping on latitude
+;
+	if (nlat eq -1) then begin
+		theloopmin=0
+		theloopmax=south_north-1
+	endif else begin
+		theloopmin=nlat-1 ;!! IDL arrays starts at 0
+		theloopmax=nlat-1
+	endelse	
+	if (when eq -1) then nt = 0 else nt=when-1
+	case vert_coord of
+		'height': charlev='Altitude (km)'
+		'pressure': begin
+			H=10.
+			p0=6.1
+			z=H*alog(p0/z)
+			charlev='Log-Pressure altitude (km)'
+		end
+		'model_level': charlev='Model level'
+	endcase
+end
+'meridional': begin
+;
+; 5. meridional: looping on longitude
+;
+	if (nlon eq -1) then begin
+		theloopmin=0
+		theloopmax=west_east-1
+	endif else begin
+		theloopmin=nlon-1 ;!! IDL arrays starts at 0
+		theloopmax=nlon-1
+	endelse	
+	if (when eq -1) then nt = 0 else nt=when-1
+	case vert_coord of
+		'height': charlev='Altitude (km)'
+		'pressure': begin
+			H=10.
+			p0=6.1
+			z=H*alog(p0/z)
+			charlev='Log-Pressure altitude (km)'
+		end
+		'model_level': charlev='Model level'
+	endcase
+end
+else: stop
+endcase
+
+
+;-------------
+; plot loop
+;-------------	
+for theloop=theloopmin,theloopmax do begin
+
+case plot_user of
+'map': begin
+	;
+	; case 1: map
+	;	
+	case vert_coord of
+	'height': charlev=string(round(z(theloop)*1000.),'(I0)')+' m'
+	'pressure': charlev=string(round(z(theloop)*100.),'(I0)')+' Pa'
+	'model_level': charlev=string(round(z(theloop)),'(I0)')
+	endcase
+	title=title_save+' at level '+charlev
+	print, title
+
+	;-------------
+	; draw map
+	;-------------
+
+	what_I_plot=var[*,*,theloop,nt]
+
+	if (diff_scale eq 0) then begin
+		minfield_init=minfield	;-same scale for all levels 
+		maxfield_init=maxfield	;-same scale for all levels
+	endif
+	if (n_elements(fieldc_user) eq 0) then begin
+		overcontour=0.
+	endif else begin
+		if (n_elements(topo) ne 0) then begin
+			overcontour=varcontour[*,*,0,0] ; other levels than 0 are blank !
+		endif else begin
+			overcontour=varcontour[*,*,theloop,nt]
+		endelse
+	endelse
+	if (n_elements(winds) eq 0) then begin
+		overvector_x=0.
+		overvector_y=0.
+	endif else begin
+		overvector_x=u[*,*,theloop,nt]
+		overvector_y=v[*,*,theloop,nt]
+
+
+;--------------------------------------------------------
+; NB: you can use the winds fields to combine variables 
+;--------------------------------------------------------
+if (n_elements(winds) eq 2) then begin
+	;** HYDROSTATIC REDUCTION TO LEVEL OF REFERENCE
+	if ((winds(0) eq 'tk') and (winds(1) eq 'HGT')) then begin
+	overvector_x=reform(u[*,*,6,nt])	;;T 1km
+	overvector_y=reform(v[*,*,0,nt])	;;topo (m)
+	zref=mean(overvector_y)
+	;zref=median(overvector_y)
+	;zref=(max(overvector_y)-min(overvector_y))/2.
+	;zref=3000.
+	print, 'zref', zref
+	H=192.*overvector_x/3.72	;unit is m
+	what_I_plot=what_I_plot*exp((overvector_y-zref)/H)
+	overvector_x=0
+	overvector_y=0
+	endif
+endif
+;--------------------------------------------------------
+; NB: you can use the winds fields to combine variables
+;--------------------------------------------------------
+
+	endelse
+
+;;******************TEMP
+;;pressure anomaly
+;nday=3
+;oneday=reform(var[*,*,theloop,0:nday-1])
+;what_I_plot=what_I_plot-(total(oneday,3)/float(nday))
+;;******************TEMP
+
+	map_latlon, $
+		what_I_plot, $
+		lon, $
+		lat, $
+		ct=pal, $
+		minfield=minfield_init, $
+		maxfield=maxfield_init, $
+		overcontour=overcontour, $
+		overvector_x=overvector_x, $
+		overvector_y=overvector_y, $
+		colors=colors, $
+		title=title, $
+		format=format	
+end
+'profile': begin
+	;
+	; case 2: profile
+	;
+	print, title
+
+	what_I_plot=var[nlon,nlat,*,theloop]
+	if (diff_scale eq 0) then begin
+		minfield_init=minfield	;-same scale for all levels 
+		maxfield_init=maxfield	;-same scale for all levels
+	endif
+	if (n_elements(fieldc_user) eq 0) then begin
+		overcontour=0.
+		overcontourz=0.
+	endif else begin
+		overcontour=varcontour[nlon,nlat,*,theloop]
+		overcontourz=0.
+	endelse
+	if (n_elements(inprofile) gt 1) then begin
+		overcontour=inprofile
+		overcontourz=0.
+		if (n_elements(incolumn) gt 1) then overcontourz=incolumn
+	endif
+
+	profile, $
+		what_I_plot, $				; 1D vertical profile
+		column, $				; altitudes	
+	;	alt=alt, $				; altitude range [altmin, altmax]
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		inprofile=overcontour, $		; another vertical profile to overplot
+		incolumn=overcontourz, $		; altitudes of the other vertical profile (in case /= column)
+		title_plot=title, $			; title of the plot ('Profile' is default)
+		title_axis=title_axis, $		; title of the [x,y] axis (['Field','Altitude'] is default)
+	;	mention=mention, $			; add text precision within the plot window (default is nothing or '')
+		discrete=discrete			; show the profile points (= type of points in !psym)
+end
+'hovmuller': begin
+	;
+	; case 3: hovmuller
+	;
+	; NB: lat/lon hovmuller 
+	;
+	what_I_plot=reform(var[*,nlat,theloop,*])
+	space=lon
+	temps=findgen(time)*2
+	if (n_elements(fieldc_user) ne 0) then begin
+		overcontour=reform(varcontour[*,nlat,theloop,*])
+	endif else begin
+		overcontour=0.
+	endelse
+	minspace=0.
+	maxspace=0.
+
+;;******************TEMP
+;	minspace=-90.
+;	maxspace=90.
+
+if (n_elements(winds) ne 0) then begin		;;;winds sert à faire des calculs !!
+overvector_x=reform(u[*,nlat,4,*])	;;T 1km
+overvector_y=reform(v[*,nlat,0,*])	;;topo (m)
+zref=median(overvector_y)
+H=192.*overvector_x/3.72	;unit is m
+;print, (overvector_y-zref)/H
+what_I_plot=what_I_plot*exp((overvector_y-zref)/H)
+
+
+endif 
+;oneday=reform(var[*,nlat,0,0:11])
+;slice=(total(oneday,2)/12.)
+;for i=0, n_elements(var[*,0,0,0])-1 do begin
+;	what_I_plot(i,*)=what_I_plot(i,*)-slice(i)
+;endfor
+;;******************TEMP
+
+
+
+
+
+	hovmuller, $
+		what_I_plot, $				; 2D field
+		space, $				; space coordinate
+		temps, $				; time coordinate
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		minspace=minspace, $			; minimum value of space window (=0: calculate)
+		maxspace=maxspace, $			; maximum value of space window (=0: calculate)
+		overcontour=overcontour, $		; another 2D field to overplot with contour lines (=0: no)
+		ct=pal, $				; color table (33-rainbow is default)
+		colors=colors, $			; number of colors/levels (32 is default)
+		title=title				; title of the plot ('' is default)
+end
+'zonal': begin
+	;
+	; case 4: zonal section
+	; 
+	displaylat=float(round(lat(theloop)*100))/100.
+	title=title_save+ $
+		' section at latitude '+string(displaylat,format='(1F6.2)')+ $
+		'!Uo!N '
+	print, title
+
+	what_I_plot=reform(var[*,theloop,*,nt])
+	space=lon
+	altitude=z
+
+	if (n_elements(fieldc_user) ne 0) then begin
+		overcontour=reform(varcontour[*,theloop,*,nt])
+	endif else begin
+		overcontour=0.
+	endelse
+
+	if (n_elements(winds) eq 0) then begin
+		overvector_x=0.
+		overvector_y=0.
+	endif else begin
+		overvector_x=reform(u[*,theloop,*,nt])		;;horizontal component
+		overvector_y=reform(v[*,theloop,*,nt])		;;vertical component
+	endelse
+
+	minspace=0.
+	maxspace=0.
+
+	title_user=title
+	title_axis=['Longitude',charlev]
+
+	section, $
+		what_I_plot, $				; 2D field
+		space, $				; horizontal coordinate
+		altitude, $				; altitude coordinate
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		minspace=minspace, $			; minimum value of space window (=0: calculate)
+		maxspace=maxspace, $			; maximum value of space window (=0: calculate)
+		overcontour=overcontour, $		; another 2D field to overplot with contour lines (=0: no)
+		overvector_x=overvector_x, $		; wind vector - x component (=0: no)
+		overvector_y=overvector_y, $		; wind vector - y component (=0: no)
+		colors=colors, $			; number of colors/levels (32 is default)
+		title_plot=title_user, $		; title of the plot ('Profile' is default)
+		title_axis=title_axis, $		; title of the [x,y] axis (['Field','Altitude'] is default)
+		ct=pal, $				; color table (33-rainbow is default)
+		format=format
+end	
+'meridional': begin
+	;
+	; case 5: meridional section
+	; 
+	displaylon=float(round(lon(theloop)*100))/100.
+	title=title_save+ $
+		' section at longitude '+string(displaylon,format='(1F7.2)')+ $
+		'!Uo!N '
+	print, title
+
+	what_I_plot=reform(var[theloop,*,*,nt])
+	space=lat
+	altitude=z
+
+	if (n_elements(fieldc_user) ne 0) then begin
+		overcontour=reform(varcontour[theloop,*,*,nt])
+	endif else begin
+		overcontour=0.
+	endelse
+
+	if (n_elements(winds) eq 0) then begin
+		overvector_x=0.
+		overvector_y=0.
+	endif else begin
+		overvector_x=reform(u[theloop,*,*,nt])		;;horizontal component
+		overvector_y=reform(v[theloop,*,*,nt])		;;vertical component
+	endelse
+
+	minspace=0.
+	maxspace=0.
+
+	title_user=title
+	title_axis=['Latitude',charlev]
+
+	section, $
+		what_I_plot, $				; 2D field
+		space, $				; horizontal coordinate
+		altitude, $				; altitude coordinate
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		minspace=minspace, $			; minimum value of space window (=0: calculate)
+		maxspace=maxspace, $			; maximum value of space window (=0: calculate)
+		overcontour=overcontour, $		; another 2D field to overplot with contour lines (=0: no)
+		overvector_x=overvector_x, $		; wind vector - x component (=0: no)
+		overvector_y=overvector_y, $		; wind vector - y component (=0: no)
+		colors=colors, $			; number of colors/levels (32 is default)
+		title_plot=title_user, $		; title of the plot ('Profile' is default)
+		title_axis=title_axis, $		; title of the [x,y] axis (['Field','Altitude'] is default)
+		ct=pal					; color table (33-rainbow is default)
+end	
+
+
+
+
+
+
+
+
+else: 
+endcase
+
+
+
+endfor
+device, /close
+
+	; save the profiles for another plot (last of the time loop)
+	if (n_elements(inprofile) le 1) and (plot_user eq 'profile') then begin
+		inprofile=what_I_plot
+		incolumn=column
+	endif
+
+;******************************************************************************
+;
+; BEGIN: different from out_wrf
+;
+;******************************************************************************
+;;;skip:
+;;;endfor
+
+end
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/out_wrf.pro
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/out_wrf.pro	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/out_wrf.pro	(revision 86)
@@ -0,0 +1,808 @@
+pro out_wrf, $
+		plot=plot_user, $		; 'map', 'profile', 'hovmuller', 'meridional', 'zonal' (default is 'map')
+	;----------------
+	; input fields
+	;----------------
+		field1=field_user, $		; field to contour: 'U', 'V', 'tk', ... etc ...
+		field2=fieldc_user, $		; field to overplot (contour lines, optional)
+	;------------------------------------------------------
+	; 'map' options (no effect when plot is 'profile')
+	;------------------------------------------------------
+		topo=topo, $			; overplot topography contours (!! must load HGT !! override field2)
+		winds=winds, $			; option to overplot winds (e.g. =['U','V'])
+		level=level, $			; vertical level subscript (starts at 1 ...)(if -1: all levels)
+	;------------------------------------------------------
+	; 'profile' options (no effect when plot is 'map')
+	;------------------------------------------------------
+		inprofile=inprofile, $		; a profile to compare with ... (override field2)
+		incolumn=incolumn, $		; altitudes ... in case /= column
+		;; NB: inprofile and incolumn may be used to output profiles
+		;; ... reset with inprofile=0, incolumn=0
+	;----------------
+	; choose time
+	;----------------
+		when=when, $			; time subscript (starts at 1 ...)(if -1: first time step)
+	;----------------
+	; choose latitude/longitude
+	;----------------
+		nlat=nlat, $			; latitude subscript (starts at 1 ...)(if -1: all latitudes)(default is middle)
+		nlon=nlon, $			; longitude subscript (starts at 1 ...)(if -1: all longitudes)(default is middle)
+	;----------------
+	; plot options
+	;----------------
+;;NB: celles ci sont dorenavant indiquees dans le _user si besoin 
+		range=range_user, $		; range=[min,max]
+		colors=colors, $		; 16, 32, 64, ... etc
+save_ps=file_ps, $
+	;-------------------------
+	; no useless data loading
+	;-------------------------
+		save_data=save_data		; to avoid loading data if already done once
+
+
+
+;---------------------------------------------------------------------------
+;
+; A general routine to plot ARWpost outputs
+;
+; USE:
+; - out_wrf	(will plot everything !)
+; - out_wrf, field1='tk', range=[160,300]
+;
+; - out_wrf, field1='pressure', field2='tk'
+; - out_wrf, field1=TSK, field2=HGT, winds=['U','V'] 
+;
+; - out_wrf, save_data=field
+;   ... then next calls including 'save_data=field' won't include the data loading 
+;
+; - out_wrf, field1='tk', winds=['umet','vmet'], topo=1, save_data=yeah, level=1, when=2
+;
+; - out_wrf, field1='tk', plot='profile'
+;
+; - out_wrf, plot='profile', field1='tk', when=1, inprofile=inprofile, incolumn=incolumn
+;   out_diagfi, plot='profile', field1='temp', when=112, inprofile=inprofile, incolumn=incolumn 	
+;
+; - out_wrf, plot='meridional',field1='U'
+;
+; - out_wrf, plot='hovmuller',field1='TSK'
+;
+; - out_wrf, plot='zonal',field1='V', nlat=-1
+; 
+;---------------------------------------------------------------------------
+; A. Spiga, August 2007, adding plot keyword September 2007
+;---------------------------------------------------------------------------
+
+
+    
+
+
+;-------------------------------------
+; default settings and user settings
+;-------------------------------------
+diff_scale=0  ;;1
+charend=''
+format=''
+if (n_elements(plot_user) eq 0) then plot_user='map'
+if (n_elements(field_user) eq 0) then field_user='all'
+if (n_elements(level) eq 0) then level=-1
+if (n_elements(when) eq 0) then when=-1
+;;if (n_elements(nlon) eq 0) then nlon=-1
+;;if (n_elements(nlat) eq 0) then nlat=-1
+if (n_elements(colors) eq 0) then colors=32
+;if (n_elements(minfield_init)*n_elements(maxfield_init) eq 0) then begin
+if (n_elements(range_user) eq 0) then begin
+	minfield_init = 0.
+	maxfield_init = 0.
+diff_scale=1
+endif
+
+if (n_elements(range_user) eq 2) then begin
+	minfield_init=float(range_user(0))
+	maxfield_init=float(range_user(1))
+	print, 'setting user limits !', minfield_init, maxfield_init 
+	diff_scale=1
+endif
+
+if (n_elements(topo) eq 0) then topo=0
+if (topo eq 1) then fieldc_user='HGT'	;case with topography
+
+if (n_elements(fieldc_user) eq 0) then fieldc_user=''
+
+;goto, squeeze
+;------------------
+; read data
+;------------------
+;
+; read the file generated by ARWpost
+;
+;*********************************************
+if (n_elements(save_data) eq 0) then begin
+	read_dat, field, x, y, z, t, desc_field
+	;save_data=field
+endif else begin
+	read_dat, field, x, y, z, t, desc_field, nodata=1  ; read only the dimensions
+	field=save_data
+	save_data=0.	;;free memory
+endelse
+;;*********************************************
+;; opt: saved IDL files
+;squeeze:
+;restore, filename='OM_nonhydro.idl'
+;field=output
+;;nh=output
+;;restore, filename='OM_hydro.idl'
+;;h=output
+;;field=abs(nh-h)
+;;*********************************************
+
+;
+; dimensions
+;
+lon=x
+west_east=n_elements(x)
+if (n_elements(nlon) eq 0) then begin
+	nlon=round(west_east/2)
+endif 
+if (nlon eq 0) then nlon=round(west_east/2)
+
+lat=y
+south_north=n_elements(y)
+if (n_elements(nlat) eq 0) then begin
+	nlat=round(south_north/2)
+endif
+if (nlat eq 0) then nlat=round(south_north/2)
+	
+bottom_top=n_elements(z)
+
+time=n_elements(t)
+if ((time eq 1) and (plot_user eq 'hovmuller')) then stop
+
+nfields=n_elements(desc_field)
+if (z(0) lt z(1)) then vert_coord='height'
+if (z(0) gt z(1)) then vert_coord='pressure'
+if ((z(0) eq 1.) and (z(bottom_top-1) eq (z(bottom_top-2)+1.))) then vert_coord='model_level'
+;print, 'vert_coord is ... ', vert_coord
+
+;
+; find the field desired by the user
+;
+if (field_user ne 'all') then begin
+	match=-1
+	indice=-1
+	repeat begin
+		indice=indice+1
+		if (indice eq nfields) $
+			then match=9999 $
+			else match=STRPOS(desc_field(indice),field_user)  
+	endrep until match ne -1
+	if (match eq 9999) then begin
+		print, field_user, ' not found !'
+		stop
+	endif
+;print, 'selected for plot:  ',desc_field(indice)
+endif
+
+;
+; prepare loop for graphics
+;
+if (n_elements(indice) ne 0) then begin
+	field_loop_start=indice
+	field_loop_end=indice
+endif else begin
+	field_loop_start=0
+	field_loop_end=nfields-1	
+endelse
+
+;
+; overplot contour of another field ?
+;
+if (fieldc_user ne '') then begin
+	match=-1
+	indice=-1
+	repeat begin
+		indice=indice+1
+		if (indice eq nfields) $
+			then match=9999 $
+			else match=STRPOS(desc_field(indice),fieldc_user)  
+	endrep until match ne -1
+	if (match eq 9999) then begin
+		print, fieldc_user, ' not found !'
+		stop
+	endif
+;print, 'selected for over-contour:  ',desc_field(indice)
+	varcontour=reform(field(indice,*,*,*,*))
+endif
+
+;
+; overplot wind vectors ?
+;
+if (n_elements(winds) ge 2) then begin
+	; zonal wind
+	match=-1
+	indice=-1
+	repeat begin
+		indice=indice+1
+		if (indice eq nfields) $
+			then match=9999 $
+			else match=STRPOS(desc_field(indice),winds(0))  
+	endrep until match ne -1
+	if (match eq 9999) then begin
+		print, 'u not found !'
+		stop
+	endif
+;print, 'selected for vector_x:  ',desc_field(indice)
+	u=reform(field(indice,*,*,*,*))
+	; meridional wind
+	match=-1
+	indice=-1
+	repeat begin
+		indice=indice+1
+		if (indice eq nfields) $
+			then match=9999 $
+			else match=STRPOS(desc_field(indice),winds(1))  
+	endrep until match ne -1
+	if (match eq 9999) then begin
+		print, 'v not found !'
+		stop
+	endif
+;print, 'selected for vector_y:  ',desc_field(indice)
+	v=reform(field(indice,*,*,*,*))
+endif
+
+
+
+;
+; loop for graphics   
+;
+for field_loop=field_loop_start,field_loop_end do begin
+
+	charstart=STRSPLIT(desc_field(field_loop), ' ', LENGTH=charlength)
+	charvar=STRMID(desc_field(field_loop),charstart(0),charlength(0))
+	var=reform(field(field_loop,*,*,*,*))
+
+	; useless field ...
+	if (charvar eq vert_coord) then goto, skip
+
+	;
+	; min/max on the whole field
+	;
+	w=where(var ne 1e37)
+	if (w(0) ne -1) then begin
+		minfield=min(var[w])
+		maxfield=max(var[w])
+			;;; manage unvalid points
+			;;;w=where(var eq 1e37)
+			;;;if (w(0) ne -1) then var[w]=!Values.F_NAN
+	endif else begin
+		var=0
+		print, 'no valid values in this field', charvar 
+	endelse
+
+;----------------------
+; initialize graphics
+;----------------------
+
+;!p.charthick = 2.0
+;!p.thick = 3.0
+;!x.thick = 2.0
+;!y.thick = 2.0
+
+
+lct = round(t[when-1]+24) MOD 24
+
+;;lct = when  ;; on choisit de les numeroter dans l'ordre 
+;;if ((when ne -1) and (lct lt 10)) then charend=charend+'0'+string(lct, '(I0)')
+;;if (lct ge 10) then charend=charend+string(lct, '(I0)') 
+;;lct = round(t[when-1]) MOD 24
+;;denom=plot_user+'_'+vert_coord+'_'+charvar
+;;if (fieldc_user ne '') then denom=denom+'_'+fieldc_user
+;;if (winds(0) ne '') then denom=denom+'_'+winds(0)+winds(1)
+;;denom=denom+charend
+
+
+if (n_elements(file_ps) ne '') then begin
+	set_plot, 'ps'
+	device, filename=file_ps+'.ps', /color
+endif
+
+
+case charvar of
+'HGT': begin
+        title='Topography (km)'
+	level=1
+	pal=16
+	format='(F7.0)'
+end
+'tk': begin
+	title='Temperature (K)'
+	pal=3
+	format='(F4.0)'
+end
+'TSK': begin
+	title='Surface Temperature (K)'
+	pal=3
+	level=1
+end
+'U': begin
+	title='Zonal wind (!Nm!N.s!U-1!N)'
+	pal=33
+	if (winds(0) eq 'U') then title='Wind field (!Nm!N.s!U-1!N)'
+end
+'V': begin
+	title='Meridional wind (!Nm!N.s!U-1!N)'
+	pal=33
+	if (winds(0) eq 'V') then title='Wind field (!Nm!N.s!U-1!N)'
+end
+'pressure': begin
+	title='Pressure (hPa)'
+	pal=16
+end
+'height': begin
+	title='Height (km)'
+	pal=16
+end
+'T': title='Potential temperature departure from t0 (K)'
+'W': begin
+	title='Vertical wind (!Nm!N.s!U-1!N)'
+	pal=16
+	format='(F6.2)'
+end
+else: begin
+	title=''
+	pal=33
+end
+endcase
+
+if (when ne -1) then begin
+	title=title+', LT '+string(lct, '(I0)')
+endif else begin
+	title=title+', time '+string(when, '(I0)')
+endelse
+	
+title_save=title
+
+
+
+;---------------------------------
+; choose positioning ...
+;---------------------------------
+
+case plot_user of
+'map': begin
+;
+; 1. map: looping on vertical levels
+;
+	if (level eq -1) then begin
+		theloopmin=0
+		theloopmax=bottom_top-1
+	endif else begin
+		theloopmin=level-1 ;!! IDL arrays starts at 0
+		theloopmax=level-1
+	endelse	
+	if (when eq -1) then nt = 0 else nt=when-1
+end
+'profile': begin
+;
+; 2. profile: looping on time
+;
+	;;discrete=4
+	title_axis=strarr(2)
+	title_axis(0)=title_save
+	if (when eq -1) then begin
+		theloopmin=0
+		theloopmax=time-1
+	endif else begin
+		theloopmin=when-1 ;!! IDL arrays starts at 0
+		theloopmax=when-1
+	endelse	
+
+	column=z
+	case vert_coord of
+		'height': title_axis(1)='Altitude (km)'
+		'pressure': begin
+			H=10.
+			p0=6.1
+			column=H*alog(p0/z)
+			title_axis(1)='Log-Pressure altitude (km)'
+		end
+		'model_level': title_axis(1)='Model level'
+	endcase
+
+	displaylat=float(round(lat(nlat)*100))/100.
+	displaylon=float(round(lon(nlon)*100))/100.
+	description='WRF Profile at latitude '+string(displaylat,format='(1F6.2)')+$
+	'!Uo'+'!N and longitude '+string(displaylon,format='(1F7.2)')+'!Uo!N '
+	title=description
+end
+'hovmuller': begin
+;
+; 3. hovmuller: looping on vertical levels
+;
+	if (level eq -1) then begin
+		theloopmin=0
+		theloopmax=bottom_top-1
+	endif else begin
+		theloopmin=level-1 ;!! IDL arrays starts at 0
+		theloopmax=level-1
+	endelse	
+end
+'zonal': begin
+;
+; 4. zonal: looping on latitude
+;
+	if (nlat eq -1) then begin
+		theloopmin=0
+		theloopmax=south_north-1
+	endif else begin
+		theloopmin=nlat-1 ;!! IDL arrays starts at 0
+		theloopmax=nlat-1
+	endelse	
+	if (when eq -1) then nt = 0 else nt=when-1
+	case vert_coord of
+		'height': charlev='Altitude (km)'
+		'pressure': begin
+			H=10.
+			p0=6.1
+			z=H*alog(p0/z)
+			charlev='Log-Pressure altitude (km)'
+		end
+		'model_level': charlev='Model level'
+	endcase
+end
+'meridional': begin
+;
+; 5. meridional: looping on longitude
+;
+	if (nlon eq -1) then begin
+		theloopmin=0
+		theloopmax=west_east-1
+	endif else begin
+		theloopmin=nlon-1 ;!! IDL arrays starts at 0
+		theloopmax=nlon-1
+	endelse	
+	if (when eq -1) then nt = 0 else nt=when-1
+	case vert_coord of
+		'height': charlev='Altitude (km)'
+		'pressure': begin
+			H=10.
+			p0=6.1
+			z=H*alog(p0/z)
+			charlev='Log-Pressure altitude (km)'
+		end
+		'model_level': charlev='Model level'
+	endcase
+end
+else: stop
+endcase
+
+
+;-------------
+; plot loop
+;-------------	
+for theloop=theloopmin,theloopmax do begin
+
+case plot_user of
+'map': begin
+	;
+	; case 1: map
+	;	
+	case vert_coord of
+	'height': charlev=string(round(z(theloop)*1000.),'(I0)')+' m'
+	'pressure': charlev=string(round(z(theloop)*100.),'(I0)')+' Pa'
+	'model_level': charlev=string(round(z(theloop)),'(I0)')
+	endcase
+	if (n_elements(var(0,0,*,0) ge 2)) then title=title_save+', at level '+charlev
+	print, title
+
+	;-------------
+	; draw map
+	;-------------
+
+	what_I_plot=var[*,*,theloop,nt]
+
+	if (diff_scale eq 0) then begin
+		minfield_init=minfield	;-same scale for all levels 
+		maxfield_init=maxfield	;-same scale for all levels
+	endif
+	if (fieldc_user eq '') then begin
+		overcontour=0.
+	endif else begin
+		if (topo eq 1) then begin
+			overcontour=varcontour[*,*,0,0]/1000.  ; other levels than 0 are blank !
+		endif else begin
+			overcontour=varcontour[*,*,theloop,nt]
+		endelse
+	endelse
+	if (n_elements(winds) lt 2) then begin
+		overvector_x=0.
+		overvector_y=0.
+	endif else begin
+		overvector_x=u[*,*,theloop,nt]
+		overvector_y=v[*,*,theloop,nt]
+
+
+;--------------------------------------------------------
+; NB: you can use the winds fields to combine variables 
+;--------------------------------------------------------
+if (n_elements(winds) ge 2) then begin
+
+
+;** HYDROSTATIC REDUCTION TO LEVEL OF REFERENCE
+@sea_level.idl
+
+;** WIND STRESS
+@wind_stress.idl
+
+;** ERTEL POTENTIAL VORTICITY
+@ertel_vorticity.idl
+
+
+endif
+;--------------------------------------------------------
+; NB: you can use the winds fields to combine variables
+;--------------------------------------------------------
+
+	endelse
+
+;;******************TEMP
+;;pressure anomaly
+;nday=3
+;oneday=reform(var[*,*,theloop,0:nday-1])
+;what_I_plot=what_I_plot-(total(oneday,3)/float(nday))
+;;******************TEMP
+
+	map_latlon, $
+		what_I_plot, $
+		lon, $
+		lat, $
+		ct=pal, $
+		minfield=minfield_init, $
+		maxfield=maxfield_init, $
+		overcontour=overcontour, $
+		overvector_x=overvector_x, $
+		overvector_y=overvector_y, $
+		colors=colors, $
+		title=title, $
+		format=format	
+end
+'profile': begin
+	;
+	; case 2: profile
+	;
+	print, title
+
+	what_I_plot=var[nlon,nlat,*,theloop]
+	if (diff_scale eq 0) then begin
+		minfield_init=minfield	;-same scale for all levels 
+		maxfield_init=maxfield	;-same scale for all levels
+	endif
+	if (fieldc_user eq '') then begin
+		overcontour=0.
+		overcontourz=0.
+	endif else begin
+		overcontour=varcontour[nlon,nlat,*,theloop]
+		overcontourz=0.
+	endelse
+	if (n_elements(inprofile) gt 1) then begin
+		overcontour=inprofile
+		overcontourz=0.
+		if (n_elements(incolumn) gt 1) then overcontourz=incolumn
+	endif
+
+	profile, $
+		what_I_plot, $				; 1D vertical profile
+		column, $				; altitudes	
+	;	alt=alt, $				; altitude range [altmin, altmax]
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		inprofile=overcontour, $		; another vertical profile to overplot
+		incolumn=overcontourz, $		; altitudes of the other vertical profile (in case /= column)
+		title_plot=title, $			; title of the plot ('Profile' is default)
+		title_axis=title_axis, $		; title of the [x,y] axis (['Field','Altitude'] is default)
+	;	mention=mention, $			; add text precision within the plot window (default is nothing or '')
+		discrete=discrete			; show the profile points (= type of points in !psym)
+end
+'hovmuller': begin
+	;
+	; case 3: hovmuller
+	;
+	; NB: lat/lon hovmuller 
+	;
+	what_I_plot=reform(var[*,nlat,theloop,*])
+	space=lon
+	temps=findgen(time)*2
+	if (fieldc_user ne '') then begin
+		overcontour=reform(varcontour[*,nlat,theloop,*])
+	endif else begin
+		overcontour=0.
+	endelse
+	minspace=0.
+	maxspace=0.
+
+;;******************TEMP
+;	minspace=-90.
+;	maxspace=90.
+
+if (n_elements(winds) ge 2) then begin		;;;winds sert à faire des calculs !!
+overvector_x=reform(u[*,nlat,4,*])	;;T 1km
+overvector_y=reform(v[*,nlat,0,*])	;;topo (m)
+zref=median(overvector_y)
+H=192.*overvector_x/3.72	;unit is m
+;print, (overvector_y-zref)/H
+what_I_plot=what_I_plot*exp((overvector_y-zref)/H)
+
+
+endif 
+;oneday=reform(var[*,nlat,0,0:11])
+;slice=(total(oneday,2)/12.)
+;for i=0, n_elements(var[*,0,0,0])-1 do begin
+;	what_I_plot(i,*)=what_I_plot(i,*)-slice(i)
+;endfor
+;;******************TEMP
+
+
+
+
+
+	hovmuller, $
+		what_I_plot, $				; 2D field
+		space, $				; space coordinate
+		temps, $				; time coordinate
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		minspace=minspace, $			; minimum value of space window (=0: calculate)
+		maxspace=maxspace, $			; maximum value of space window (=0: calculate)
+		overcontour=overcontour, $		; another 2D field to overplot with contour lines (=0: no)
+		ct=pal, $				; color table (33-rainbow is default)
+		colors=colors, $			; number of colors/levels (32 is default)
+		title=title				; title of the plot ('' is default)
+end
+'zonal': begin
+	;
+	; case 4: zonal section
+	; 
+	displaylat=float(round(lat(theloop)*100))/100.
+	title=title_save+ $
+		' section at latitude '+string(displaylat,format='(1F6.2)')+ $
+		'!Uo!N '
+	print, title
+
+	what_I_plot=reform(var[*,theloop,*,nt])
+	space=lon
+	altitude=z
+
+if (topo eq 1) then begin
+	topography=varcontour[theloop,*,0,0] ; other levels than 0 are blank !
+        fieldc_user=''
+endif else begin
+	topography=0.
+endelse	
+		    
+	
+	if (fieldc_user ne '') then begin
+		overcontour=reform(varcontour[*,theloop,*,nt])
+	endif else begin
+		overcontour=0.
+	endelse
+
+	if (n_elements(winds) lt 2) then begin
+		overvector_x=0.
+		overvector_y=0.
+	endif else begin
+		overvector_x=reform(u[*,theloop,*,nt])		;;horizontal component
+		overvector_y=reform(v[*,theloop,*,nt])		;;vertical component
+	endelse
+
+	minspace=0.
+	maxspace=0.
+
+	title_user=title
+	title_axis=['Longitude',charlev]
+
+	section, $
+		what_I_plot, $				; 2D field
+		space, $				; horizontal coordinate
+		altitude, $				; altitude coordinate
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		minspace=minspace, $			; minimum value of space window (=0: calculate)
+		maxspace=maxspace, $			; maximum value of space window (=0: calculate)
+		overcontour=overcontour, $		; another 2D field to overplot with contour lines (=0: no)
+		overvector_x=overvector_x, $		; wind vector - x component (=0: no)
+		overvector_y=overvector_y, $		; wind vector - y component (=0: no)
+		colors=colors, $			; number of colors/levels (32 is default)
+		title_plot=title_user, $		; title of the plot ('Profile' is default)
+		title_axis=title_axis, $		; title of the [x,y] axis (['Field','Altitude'] is default)
+		ct=pal, $				; color table (33-rainbow is default)
+topo=topography, $
+		format=format
+end	
+'meridional': begin
+	;
+	; case 5: meridional section
+	; 
+	displaylon=float(round(lon(theloop)*100))/100.
+	title=title_save+ $
+		' section at longitude '+string(displaylon,format='(1F7.2)')+ $
+		'!Uo!N '
+	print, title
+
+	what_I_plot=reform(var[theloop,*,*,nt])
+	space=lat
+	altitude=z
+
+if (topo eq 1) then begin
+        topography=varcontour[theloop,*,0,0] ; other levels than 0 are blank !
+        fieldc_user=''
+endif else begin
+        topography=0.
+endelse
+			
+	if (fieldc_user ne '') then begin
+		overcontour=reform(varcontour[theloop,*,*,nt])
+	endif else begin
+		overcontour=0.
+	endelse
+
+	if (n_elements(winds) lt 2) then begin
+		overvector_x=0.
+		overvector_y=0.
+	endif else begin
+		overvector_x=reform(u[theloop,*,*,nt])		;;horizontal component
+		overvector_y=reform(v[theloop,*,*,nt])		;;vertical component
+	endelse
+
+	minspace=0.
+	maxspace=0.
+
+	title_user=title
+	title_axis=['Latitude',charlev]
+
+	section, $
+		what_I_plot, $				; 2D field
+		space, $				; horizontal coordinate
+		altitude, $				; altitude coordinate
+		minfield=minfield_init, $		; minimum value of plotted field (=0: calculate)
+		maxfield=maxfield_init, $		; maximum value of plotted field (=0: calculate)
+		minspace=minspace, $			; minimum value of space window (=0: calculate)
+		maxspace=maxspace, $			; maximum value of space window (=0: calculate)
+		overcontour=overcontour, $		; another 2D field to overplot with contour lines (=0: no)
+		overvector_x=overvector_x, $		; wind vector - x component (=0: no)
+		overvector_y=overvector_y, $		; wind vector - y component (=0: no)
+		colors=colors, $			; number of colors/levels (32 is default)
+		title_plot=title_user, $		; title of the plot ('Profile' is default)
+		title_axis=title_axis, $		; title of the [x,y] axis (['Field','Altitude'] is default)
+topo=topography, $
+		ct=pal, $				; color table (33-rainbow is default)
+		format=format
+end	
+
+
+
+
+
+
+
+
+else: 
+endcase
+
+
+print, '...ok'
+endfor
+if (n_elements(file_ps) ne '') then device, /close
+
+	; save the profiles for another plot (last of the time loop)
+	if (n_elements(inprofile) le 1) and (plot_user eq 'profile') then begin
+		inprofile=what_I_plot
+		incolumn=column
+	endif
+
+skip:
+endfor
+
+save_data=field
+field=0.	;; free memory
+
+end
+
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/plot.def
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/plot.def	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/plot.def	(revision 86)
@@ -0,0 +1,159 @@
+;;---------------------------------------
+;
+;  fill here plot parameters then type:
+;
+;  @trace.idl
+;  -or- @report.idl
+;  -or- @movie.idl
+;
+;  this file will be saved !
+;
+;;---------------------------------------
+
+
+;-------------------------------
+;-------------------------------
+; namelist.ARWPOST parameters
+;-------------------------------
+;-------------------------------
+
+;; DEFAULT: ARWPOST.EXE ALREADY RUN (i.e. input.ctl and input.dat already there)
+;; UNCOMMENT FOLLOWING LINE TO EXTRACT FILES 
+
+;extract='yes'	;;   EXECUTE ARWPOST.EXE
+		;;   WITH THE FOLLOWING PARAMETERS :
+
+nam1=" fields = 'tk,HGT,PSFC,U,V,W' "		;; CHANGE HERE FIELDS TO BE EXTRACTED 
+nam1=" fields = 'tk' "
+nam1=" fields = 'HGT,U,V' "
+nam1=" fields = 'V,W' "
+nam1=" fields = 'HGT,ICETOT,RAVE,QH2O_ICE' "
+
+nam2=" start_date = '2024-04-04_06:00:00' "	;; CHANGE HERE START DATE
+nam3=" end_date   = '2024-04-04_08:00:00' "	;; CHANGE HERE END DATE
+nam4=" interval_seconds = 14800 "		;; SET INTERVAL   2h: 7400s, 4h: 14800s, 6h: 22200s	
+nam5=" interp_method = 1 "			;; SET VERT. COORD. 0: eta levels, 1: height (inc) or pressure (dec)
+tabnam=[1,30,0.5]   	;; SET HEIGHT/PRESSURE LEVELS WITH [first alt, last alt, step]
+
+
+
+;----------------------------------------------------
+;----------------------------------------------------
+; what do you want to plot ?
+; THIS MUST BE SET ACCORDINGLY TO PREVIOUS SECTION
+;   - do not plot a field not extracted
+;   - do not ask for more plot than data available
+;----------------------------------------------------
+;----------------------------------------------------
+
+num=1		   	;; how many plots ?
+num=6                   ;; how many plots ?
+;num=2                   ;; how many plots ?
+num=3
+
+;@map_sea_level.idl
+;@map_ertel_vorticity.idl
+;@map_wind_stress.idl
+;field1='HGT'
+;field2='U'
+;winds=['U','V']
+;topo=1
+
+;winds=['U','V'] & field1='U' & topo=1
+;winds=['V','W'] & field1='V' 
+field1='ICETOT' & topo=1
+;field1='QH2O_ICE'
+
+
+
+;----------------------------------
+;----------------------------------
+; what kind of plot do you want ? 
+;----------------------------------
+;----------------------------------
+
+plot='map' & level=-1	
+level=2		; first level is 1
+level=24
+;level=10
+;level=1
+
+;plot='meridional' & nlon=-1 	
+;nlon=75
+
+;plot='zonal' & nlat=-1		
+;nlat=100
+
+
+;----------------------------
+;----------------------------
+; change plot settings 
+;----------------------------
+;----------------------------
+
+;; uncomment to modify default plot settings
+;; | | | | | | | | | | | | | | | | | | | | |
+;; v v v v v v v v v v v v v v v v v v v v v
+
+minfield_init=0.		;; grepALL
+maxfield_init=16.		;; grepALL
+;intervalx=5.	           	;; grepALL
+;intervaly=5.	           	;; grepALL
+;title_user='my own plot'	;; grepALL
+;title_axis=['abs','ord']	;; grepALL
+;colors=30			;; grepALL
+;pal=33				;; grepALL
+;pal=0				;; grepALL
+;pal=22				;; grepALL
+;format='(I0)'			;; grepALL
+;windex=20.			;; grepALL
+;stride=2			;; grepALL
+;flag_cb='false'		;; grepALL
+;what_I_plot=what_I_plot*1e6    ;; grepALL
+
+;minalt=-4.			;; grepSEC
+;maxalt=1.5			;; grepSEC
+;;altitude=altitude*1000.  	;; grepSEC
+;factor=30.			   ;; grepSEC
+;space=space*58. & spacekm='true'   ;; grepSEC
+;minspace=10.4*58./factor             ;; grepSEC
+;maxspace=12.*58./factor            ;; grepSEC
+
+windowx=[-145,-100]	;; grepMAP
+windowy=[-25,25]	;; grepMAP
+;lon=lon*58.		;; grepMAP
+;lat=lat*58.		;; grepMAP
+
+;lim_max = 0.03 				;; 1/3 grepALL
+;w=where(what_I_plot ge lim_max) 		;; 2/3 grepALL
+;if (w[0] ne -1) then what_I_plot[w]=lim_max	;; 3/3 grepALL
+
+;lim_blank = 10 					;; 1/3 grepALL
+;w=where(abs(what_I_plot) le lim_blank)			;; 2/3 grepALL
+;if (w[0] ne -1) then what_I_plot[w]=missing_value	;; 3/3 grepALL
+
+;title_user='Sea-level surface pressure'+title_user	;; grepALL
+
+;what_I_plot=what_I_plot*1000.  				;; grepALL
+;title_user='Water ice column (!7l!3m-pr)'+title_user	;; grepALL
+subtitle_user='LMD mesoscale model'     		;; grepALL
+
+;what_I_plot(*,*)=missing_value		;; grepMAP
+
+
+;------------------------------
+;------------------------------
+; some last possible settings 
+;------------------------------
+;------------------------------
+
+backup_data='no'   	;; save ARWpost-generated text files ? default is no
+already_data='no'  	;; load ARWpost-generated text files ? default is no
+datafolder='./WHERE_RESULTS_ARE/'	;; default is ./
+plotfolder='./WHERE_PLOTS_ARE/'	        ;; default is ./
+
+textitle='TH20_WC'	;; one word.  default is 'plot'
+texcomments='sol4'	;; one word.  default is 'plot'
+
+
+
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/read_dat.pro
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/read_dat.pro	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/read_dat.pro	(revision 86)
@@ -0,0 +1,220 @@
+pro read_dat, output, x, y, z, t, desc_field, nodata=nodata
+
+
+;---------------------------------------------------------------------------
+;
+; Read the files generated by ARWpost
+;
+; Please create symbolic link 'input.ctl' and 'input.dat'
+; --- or rename files
+;
+; - output contains all fields
+; - x,y,z,t contains coordinates
+; - desc_field contains strings describing fields
+; - if nodata is not blank, do not read the fields
+;
+;---------------------------------------------------------------------------
+; A. Spiga, August 2007
+;---------------------------------------------------------------------------
+
+
+;
+; info for localtime
+;
+interv=0.
+openr, 99, 'timefil'
+readf,99,interv
+close, 99
+
+
+
+; get preliminar info on coordinates
+;
+wnx=1
+wny=1
+wnz=1
+wnt=1
+txt=''
+param=[0.,0.,0.,0.,0.,0.,0.,0.]
+openr, 1, 'input.ctl'
+while not eof(1) do begin
+    readf,1,txt
+    match=STRPOS(txt,'linear')
+	if (match ne -1) then begin
+		wnx=(STRPOS(txt,'xdef') eq -1)*wnx  ;0 if linear
+		wny=(STRPOS(txt,'ydef') eq -1)*wny  ;0 if linear
+                wnz=(STRPOS(txt,'zdef') eq -1)*wnz  ;0 if linear 
+                wnt=(STRPOS(txt,'tdef') eq -1)*wnt  ;0 if linear
+                paramstart=STRSPLIT(txt, 'r', LENGTH=paramlength)
+        	paramchar=STRMID(txt,paramstart(1),paramlength(1))
+		READS, paramchar, start, step
+		param=param+[start*(STRPOS(txt,'xdef') ne -1),$
+			step*(STRPOS(txt,'xdef') ne -1),$
+			start*(STRPOS(txt,'ydef') ne -1),$
+			step*(STRPOS(txt,'ydef') ne -1),$
+			start*(STRPOS(txt,'zdef') ne -1),$
+			step*(STRPOS(txt,'zdef') ne -1),$
+			start*(STRPOS(txt,'tdef') ne -1),$
+			step*(STRPOS(txt,'tdef') ne -1)]
+	endif	
+endwhile
+close, 1
+timebegin=param(6)
+	
+;---------------------------------
+; READ .CTL FILE INFORMATIONS
+;---------------------------------
+
+info=read_ascii('input.ctl', missing_value=1e37)
+infodat=info.field1
+
+;
+; second column : dimensions
+;
+
+w=where((infodat(0,*) eq 1e37) and (infodat(1,*) lt 1e37))
+dimensions=reform(infodat(1,w))
+
+nx=round(dimensions(0)) 
+ny=round(dimensions(1)) 
+nz=round(dimensions(2)) 
+nt=round(dimensions(3)) 
+fields=round(dimensions(4)) 
+vertical=intarr(fields)
+for i=0,fields-1 do begin
+	vertical(i)=round(dimensions(5+i))
+endfor
+
+
+;
+; first column : coordinates
+;
+
+w=where((infodat(1,*) eq 1e37) and (infodat(0,*) lt 1e37))
+if (w(0) ne -1) then coordinates=reform(infodat(0,w))
+
+if (wnx ne 0) then begin
+	wnx=nx
+	x=coordinates(0:nx-1)
+endif else begin
+	xmin=param(0)
+	xmax=param(0)+param(1)*(nx-1) MOD 360
+	;inc=param(0)+param(1)*(nx-1)	;cas mercator
+	;x=start+inc*findgen(nx) 
+	x=xmin+(xmax-xmin)*findgen(nx)/(nx-1)
+	;print, x
+	;stop
+endelse
+if (wny ne 0) then begin
+        wny=ny
+	y=coordinates(wnx:wnx+wny-1)
+endif else begin
+        start=param(2)
+        step=param(3)
+        y=start+step*findgen(ny)
+endelse
+if (wnz ne 0) then begin
+        wnz=nz
+	z=coordinates(wnx+wny:wnx+wny+wnz-1)
+endif else begin
+        start=param(4)
+        step=param(5)
+        z=start+step*findgen(nz)
+endelse
+t=interv*findgen(nt)/3700. + timebegin + mean(x)/15.
+	
+
+;
+; get info on fields
+;
+openr, 1, 'input.ctl'
+toto=5+wnx+1+wny+1+wnz+2+fields
+infofields=strarr(toto)
+readf, 1, infofields
+close, 1
+
+;
+desc_field=transpose(infofields(toto-fields:toto-1))
+
+
+if (n_elements(nodata) eq 0) then begin
+;---------------------------------
+; READ .DAT FILE DATA
+;---------------------------------
+; GrADS views gridded data sets as 5-dimensional arrays 
+; varying in longitude, latitude, vertical level, variable, and time.
+;---------------------------------
+
+print, '----------------'
+print, 'reading data ...'
+print, '----------------'
+print, 'longitudes: ', min(x), max(x), nx, ' points'
+print, 'latitudes: ', min(y), max(y), ny, ' points' 
+print, 'altitudes: ', min(z), max(z), nz, ' points'
+print, 'local times: ', min(t), max(t), nt, ' points'
+print, fields, ' recorded fields :'
+print, desc_field
+print, '----------------'
+
+	;;very slow method
+	;;----------------	
+	;data=read_ascii('input.dat', missing_value=1e37)
+	;datadat=data.field1
+
+; see http://www.dfanning.com/tips/ascii_column_data.html	
+OPENR, lun, 'input.dat', /GET_LUN
+
+;count=0L
+output=fltarr(fields,nx,ny,nz,nt)
+dummy=0.
+
+for l=0,nt-1 do begin 						;time loop
+print, 'time loop ...', l+1, ' / ', string(nt,'(I0)')  	
+	for f=0, fields-1 do begin ;fields (whether 2D or 3D)	;variable loop
+	nzf=vertical(f)
+	print, 'reading field ...', f+1, ' / ', string(fields,'(I0)') 
+		for k=0,nzf-1 do begin  			;vertical loop
+			for j=0,ny-1 do begin			;latitude loop
+			for i=0,nx-1 do begin			;longitude loop
+
+				;;very slow method
+				;output(f,i,j,k,l)=datadat(count)
+                                ;count=count+1
+
+				READF, lun, dummy
+				output(f,i,j,k,l)=dummy
+
+			endfor
+			endfor
+		endfor
+	endfor
+endfor
+
+; missing values
+w=where(output eq 1e20)
+if (w(0) ne -1) then begin
+	print, 'change missing values ...'
+	output[w]=1e37
+endif
+
+
+CLOSE, lun
+;SPAWN, '\rm input.dat_tmp'
+print, '...done !'
+help, /memory
+print, '---------'
+
+;if (count ne n_elements(datadat)) then begin
+;	print, 'some data were not read !'
+;	stop
+;endif
+
+;;
+;; Manage missing data so as IDL recognize them
+;;
+;w=where(output eq 1e37)
+;if (w(0) ne -1) then output[w]=!Values.F_NAN
+
+endif
+
+end
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/report.idl
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/report.idl	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/report.idl	(revision 86)
@@ -0,0 +1,87 @@
+;;---------------------------------------
+;  REPORT
+;
+;  **** DO NOT MODIFY ****
+;
+;  set parameters in plot.def then type:
+;
+;  @report.idl
+;
+;;---------------------------------------
+
+
+field1='' & field2='' & winds=''
+backup_data='no' & already_data='no'  	
+datafolder='./' & plotfolder='./'
+topo=0
+extract='no'
+
+textitle='plot'
+texcomments='plot'	
+
+SPAWN, '\rm user.idl'
+SPAWN, 'grep -v grep plot.def > user.idl'
+
+;; to avoid dumb and useless warnings from IDL
+space=0. & altitude=0. & lon=0. & lat=0.
+what_I_plot=0. & title_user='' & missing_value=0.
+
+@user.idl
+if (nam5 eq " interp_method = 0 ") then coord='model_level'
+if ((nam5 eq " interp_method = 1 ") and (tabnam(0) lt tabnam(1))) then coord='height'
+
+.compile call_arwpost
+if (extract eq 'yes') then call_arwpost, nam1, nam2, nam3, nam4, nam5, tabnam
+
+denom=plot+'_'+coord+'_'+field1
+if ((field2 ne '') and (topo eq 0)) then denom=denom+'_'+field2
+if (topo eq 1) then denom=denom+'_HGT'
+if (winds(0) ne '') then denom=denom+'_'+winds(0)+winds(1)
+
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.ctl input.ctl'
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.dat input.dat'
+
+SPAWN, '\rm param_plot.idl'
+SPAWN, 'grep grepALL plot.def > param_plot.idl'
+if ((plot eq 'meridional') or (plot eq 'zonal')) then SPAWN, 'grep grepSEC plot.def >> param_plot.idl'
+if (plot eq 'map') then SPAWN, 'grep grepMAP plot.def >> param_plot.idl'
+.compile map_latlon
+.compile section
+
+!p.charthick = 2.0
+!p.thick = 3.0
+!x.thick = 2.0
+!y.thick = 2.0
+
+for i=1,num do $
+        out_wrf, $
+        plot=plot, $
+        field1=field1, $
+        field2=field2, $
+        when=i, $
+        winds=winds, $
+        level=level,$
+        save_data=yeah,$        ;; data is loaded only once
+        topo=topo,$
+;        range=range,$
+        colors=colors, $
+        nlon=nlon, $
+        nlat=nlat, $
+        save_ps=denom+string(1000+i,'(I0)')
+
+SPAWN, 'touch report.comments'
+SPAWN, '\rm report.comments'
+SPAWN, 'echo '+textitle+' > report.comments'
+SPAWN, 'echo '+texcomments+' >> report.comments'
+SPAWN, 'rm -f '+denom+'_.ps'
+print, 'executing ... ','base '+denom+'????.ps < report.comments'
+SPAWN, 'base '+denom+'????.ps < report.comments'
+SPAWN, 'mv -f base.ps '+denom+'_.ps'
+
+SPAWN, 'mv -f ./'+denom+'????.ps '+plotfolder+'/'
+SPAWN, 'mv -f ./'+denom+'_.ps '+plotfolder+'/'
+if (backup_data eq 'yes') then SPAWN, 'cp -f namelist.ARWpost '+datafolder+'/'+coord+'.namelist.ARWpost'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.ctl '+datafolder+'/'+coord+'.input.ctl'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.dat '+datafolder+'/'+coord+'.input.dat'
+
+SPAWN, 'cp -f plot.def '+plotfolder+'/'+denom+'.def'
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/see.idl
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/see.idl	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/see.idl	(revision 86)
@@ -0,0 +1,76 @@
+;;---------------------------------------
+;  TRACE
+;
+;  **** DO NOT MODIFY ****
+;
+;  set parameters in plot.def then type:
+;
+;  @trace.idl
+;
+;;---------------------------------------
+
+field1='' & field2='' & winds=''
+backup_data='no' & already_data='no'  	
+datafolder='./' & plotfolder='./'
+extract='no'
+
+SPAWN, '\rm user.idl'
+SPAWN, 'grep -v grep plot.def > user.idl'
+
+;; to avoid dumb and useless warnings from IDL
+space=0. & altitude=0. & lon=0. & lat=0.  
+what_I_plot=0. & title_user='' & missing_value=0.
+
+@user.idl
+if (nam5 eq " interp_method = 0 ") then coord='model_level'
+if ((nam5 eq " interp_method = 1 ") and (tabnam(0) lt tabnam(1))) then coord='height'
+
+.compile call_arwpost
+if (extract eq 'yes') then call_arwpost, nam1, nam2, nam3, nam4, nam5, tabnam
+
+denom=plot+'_'+coord+'_'+field1
+if ((field2 ne '') and (topo eq 0)) then denom=denom+'_'+field2
+if (topo eq 1) then denom=denom+'_HGT'
+if (winds(0) ne '') then denom=denom+'_'+winds(0)+winds(1)
+
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.ctl input.ctl'
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.dat input.dat'
+
+SPAWN, '\rm param_plot.idl'
+SPAWN, 'grep grepALL plot.def > param_plot.idl'
+if ((plot eq 'meridional') or (plot eq 'zonal')) then SPAWN, 'grep grepSEC plot.def >> param_plot.idl'
+if (plot eq 'map') then SPAWN, 'grep grepMAP plot.def >> param_plot.idl'
+.compile map_latlon
+.compile section
+
+
+set_plot, 'x'
+!P.MULTI = [0, 2, 3]
+window, 0, xsize=600, ysize=900
+
+!p.charthick = 2.0
+!p.thick = 3.0
+!x.thick = 2.0
+!y.thick = 2.0
+
+
+
+for i=1,num do $
+	out_wrf, $
+	plot=plot, $
+	field1=field1, $
+	field2=field2, $
+	when=i, $
+	winds=winds, $
+	level=level,$
+        save_data=yeah,$	;; data is loaded only once
+	topo=topo,$
+;	range=range,$
+	colors=colors, $
+	nlon=nlon, $
+	nlat=nlat
+
+if (backup_data eq 'yes') then SPAWN, 'cp -f namelist.ARWpost '+datafolder+'/'+coord+'.namelist.ARWpost'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.ctl '+datafolder+'/'+coord+'.input.ctl'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.dat '+datafolder+'/'+coord+'.input.dat'
+SPAWN, 'cp -f plot.def '+plotfolder+'/'+denom+'.def'
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/trace.idl
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/trace.idl	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/trace.idl	(revision 86)
@@ -0,0 +1,93 @@
+;;---------------------------------------
+;  TRACE
+;
+;  **** DO NOT MODIFY ****
+;
+;  set parameters in plot.def then type:
+;
+;  @trace.idl
+;
+;;---------------------------------------
+
+field1='' & field2='' & winds=''
+backup_data='no' & already_data='no'  	
+datafolder='./' & plotfolder='./'
+topo=0
+extract='no'
+
+SPAWN, '\rm user.idl'
+SPAWN, 'grep -v grep plot.def > user.idl'
+
+;; to avoid dumb and useless warnings from IDL
+space=0. & altitude=0. & lon=0. & lat=0.  
+what_I_plot=0. & title_user='' & missing_value=0.
+
+@user.idl
+if (nam5 eq " interp_method = 0 ") then coord='model_level'
+if ((nam5 eq " interp_method = 1 ") and (tabnam(0) lt tabnam(1))) then coord='height'
+
+.compile call_arwpost
+if (extract eq 'yes') then call_arwpost, nam1, nam2, nam3, nam4, nam5, tabnam
+
+denom=plot+'_'+coord+'_'+field1
+if ((field2 ne '') and (topo eq 0)) then denom=denom+'_'+field2
+if (topo eq 1) then denom=denom+'_HGT'
+if (winds(0) ne '') then denom=denom+'_'+winds(0)+winds(1)
+
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.ctl input.ctl'
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+datafolder+'/'+coord+'.input.dat input.dat'
+
+SPAWN, '\rm param_plot.idl'
+SPAWN, 'grep grepALL plot.def > param_plot.idl'
+if ((plot eq 'meridional') or (plot eq 'zonal')) then SPAWN, 'grep grepSEC plot.def >> param_plot.idl'
+if (plot eq 'map') then SPAWN, 'grep grepMAP plot.def >> param_plot.idl'
+.compile map_latlon
+.compile section
+
+!p.charthick = 2.0
+!p.thick = 3.0
+!x.thick = 2.0
+!y.thick = 2.0
+
+out_wrf, $
+        plot=plot, $
+        field1=field1, $
+        field2=field2, $
+        when=1, $
+        winds=winds, $
+        level=level,$
+        save_data=yeah,$        ;; data is loaded only once
+        topo=topo,$
+;        range=range,$
+        colors=colors, $
+        nlon=nlon, $
+        nlat=nlat, $
+        save_ps=denom+string(1001,'(I0)')
+;
+; display the first one as a sample
+;
+SPAWN, 'display ./'+denom+'1001.ps &'
+
+if (num ge 2) then for i=2,num do $
+	out_wrf, $
+	plot=plot, $
+	field1=field1, $
+	field2=field2, $
+	when=i, $
+	winds=winds, $
+	level=level,$
+        save_data=yeah,$	;; data is loaded only once
+	topo=topo,$
+	range=range,$
+	colors=colors, $
+	nlon=nlon, $
+	nlat=nlat, $
+	save_ps=denom+string(1000+i,'(I0)')
+
+SPAWN, 'mv -f ./'+denom+'????.ps '+plotfolder+'/'
+if (backup_data eq 'yes') then SPAWN, 'cp -f namelist.ARWpost '+datafolder+'/'+coord+'.namelist.ARWpost'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.ctl '+datafolder+'/'+coord+'.input.ctl'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.dat '+datafolder+'/'+coord+'.input.dat'
+
+SPAWN, 'cp -f plot.def '+plotfolder+'/'+denom+'.def'
+
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/trace.idl_old
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/trace.idl_old	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/trace.idl_old	(revision 86)
@@ -0,0 +1,45 @@
+;;---------------------------------------
+;  TRACE
+;
+;  set parameters in user.idl then type:
+;
+;  @trace.idl
+;
+;;---------------------------------------
+
+@user.idl
+
+denom=plot+'_'+coord+'_'+field1
+
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+folder+'/'+coord+'.input.ctl input.ctl'
+if (already_data eq 'yes') then SPAWN, 'ln -sf  '+folder+'/'+coord+'.input.dat input.dat'
+
+SPAWN, '\rm map_latlon_user'
+SPAWN, 'tail -15 user.idl > map_latlon_user'
+.compile map_latlon
+
+SPAWN, '\rm section_user'
+SPAWN, 'tail -15 user.idl > section_user'
+.compile section
+
+for i=1,num do $
+	out_wrf, $
+	plot=plot, $
+	field1=field1, $
+	field2=field2, $
+	when=i, $
+	winds=winds, $
+	level=level,$
+        save_data=yeah,$
+	topo=topo,$
+	range=range,$
+	colors=colors, $
+	nlon=nlon, $
+	nlat=nlat 
+
+SPAWN, 'mv -f ./'+denom+'??.ps '+folder+'/'
+if (backup_data eq 'yes') then SPAWN, 'cp -f namelist.ARWpost '+folder+'/'+coord+'.namelist.ARWpost'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.ctl '+folder+'/'+coord+'.input.ctl'
+if (backup_data eq 'yes') then SPAWN, 'cp -f input.dat '+folder+'/'+coord+'.input.dat'
+
+SPAWN, 'cp -f user.idl '+folder+'/'+denom+'.idl'
Index: trunk/mesoscale/PLOT/RESERVE/obsolete/user.idl_old
===================================================================
--- trunk/mesoscale/PLOT/RESERVE/obsolete/user.idl_old	(revision 86)
+++ trunk/mesoscale/PLOT/RESERVE/obsolete/user.idl_old	(revision 86)
@@ -0,0 +1,126 @@
+;;---------------------------------------
+;
+;  fill here plot parameters then type:
+;
+;  @trace.idl
+;  @report.idl
+;  @movie.idl
+;
+;;---------------------------------------
+
+;-----------------
+; initialization
+;-----------------
+winds=0.
+range=0.
+topo=0.
+colors=32
+backup_data='no' ;backup_data='yes'
+already_data='no' ;already_data='yes'
+folder='./WHERE_RESULTS_ARE/'
+
+
+;-----------------
+; check this first
+;-----------------
+num=6
+;
+;coord='height'
+coord='model_level'
+
+;----------------------------
+; general plot settings
+;----------------------------
+plot='map' & level=-1	
+plot='meridional' & nlon=-1 	;bug a corriger: 0
+;plot='zonal' & nlat=-1		;bug a corriger: 0
+
+;-------------------------------------------
+; easy access to already defined plots
+;-------------------------------------------
+;@map_sea_level.idl
+@map_winds_only.idl
+;@map_ertel_vorticity.idl
+;@map_wind_stress.idl
+
+;------------------------
+; user-defined settings
+;------------------------
+; keywords : plot, winds, field1, field2, range, topo, colors
+;------------------------
+
+;field1='tk'
+;nlon=100
+;range=[200,250]
+
+
+
+
+title_user=''     ;to avoid a warning later
+what_I_plot=0.    ;to avoid a warning later
+missing_value=0.  ;to avoid a warning later
+space=0.
+;;***************************************************
+;; below this line do not modify number of lines
+;; because this will be copied ... with command tail
+;; comment a line if you do not need it :)
+;;***************************************************
+;;------------------------
+;; modify plot parameters
+
+;; SECTION: 
+title_user, title_axis (2) 
+        xtitle='Horizontal coordinate'
+        ytitle='Altitude'
+endif else begin
+        xtitle=title_axis(0)
+        ytitle=title_axis(1)
+endelse
+
+if (n_elements(space) eq 0) then begin
+        space=findgen(n_elements(what_I_plot(*,0)))
+endif
+if (n_elements(altitude) eq 0) then begin
+        altitude=findgen(n_elements(what_I_plot(0,*)))
+endif
+
+if (n_elements(minfield_init) eq 0) then minfield_init=0.
+if (n_elements(maxfield_init) eq 0) then maxfield_init=0.
+if (n_elements(minspace) eq 0) then minspace=0.
+if (n_elements(maxspace) eq 0) then maxspace=0.
+if (n_elements(minalt) eq 0) then minalt=min(altitude)
+if (n_elements(maxalt) eq 0) then maxalt=max(altitude)
+if (n_elements(colors) eq 0) then colors=32
+if (n_elements(title) eq 0) then title=''
+if (n_elements(pal) eq 0) then pal=33
+if (n_elements(format) eq 0) then format='(F6.2)'
+if (format eq '') then format='(F6.2)'
+
+if ((n_elements(overvector_x) eq 0) or (n_elements(overvector_y) eq 0)) then begin
+        overvector_x=0.
+        overvector_y=0.
+endif
+
+;
+; dilatation factor for winds
+;
+if (n_elements(overvector_x)*n_elements(overvector_y) ne 1) then space=space/factor
+
+
+
+
+
+
+;;------------------------
+;;
+space=space 
+ ;minspace=158.5 & maxspace=160.
+minalt=-3 & maxalt=7
+ ;factor=15.
+ ;windex=20.
+ ;format='(F6.3)'
+ ;title_user='Sea-level surface pressure'+title_user
+pal=33
+ ;w=where(abs(what_I_plot) ge 0.03)
+ ;if (w[0] ne -1) then what_I_plot[w]=0.03
+;;**** END ****
Index: trunk/mesoscale/PLOT/SPEC/LES/getget.pro
===================================================================
--- trunk/mesoscale/PLOT/SPEC/LES/getget.pro	(revision 85)
+++ 	(revision )
@@ -1,125 +1,0 @@
-function getget, name, charvar, count=count, offset=offset, stride=stride, coordmean=coordmean, anomaly=anomaly
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-; name: char
-; charvar: char
-; count: 4 elements array   ; The dimensions of the sub-image
-; offset: 4 elements array  ; A variable that contains the offset for the sub-image
-; stride: 4 elements array  ; Create a variable to be used as a value for the STRIDE keyword
-; coordmean: 1 to 4 elements array (dimension 1 is 1)
-; anomaly: a non-null element that would be filled with mean
-;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
-
-;print, charvar
-
-
-nosave = 1  ;; comment to get an IDL data file
-notime = 1  ;; comment to get mention of time taken
-if (n_elements(notime) eq 0) then timetime = SYSTIME(1)
-
-;
-; MINIMAL ARGUMENTS
-;
-if (n_elements(name) eq 0) then stop
-if (n_elements(charvar) eq 0) then stop
-
-;
-; IDL DATA DIRECTORY
-;
-name_idl = charvar+name+'.idl'
-  if (n_elements(coordmean) ne 0) then begin
-    nc = n_elements(coordmean)
-    for i=0,nc-1 do begin
-     name_idl = string(coordmean(i),'(I0)')+name_idl   
-    endfor
-  endif
-name_idl = './idl_dat/'+name_idl
-
-;
-; TEST IF FILE EXISTS
-;
-openr,unit,name_idl,/get_lun,error=err
-
-;
-; IF IDL FILE EXISTS, THAT IS RATHER SIMPLE
-;
-IF (err EQ 0) THEN BEGIN
-
-   print, 'read from file'
-   free_lun,unit
-   restore, filename=name_idl
-
-;
-; IF IDL FILE DOES NOT EXIST, WELL, READ NETCDF
-;
-ENDIF ELSE BEGIN
-
-  ;
-  ; check save directory
-  ;
-  if (n_elements(nosave) eq 0) then SPAWN, 'mkdir -p idl_dat'  
-
-  ;
-  ; open file
-  ;
-  id=ncdf_open(name) 
-
-  ;
-  ; get dimensions
-  ;
-  var = intarr(4)
-  NCDF_DIMINQ, id, NCDF_DIMID(id, 'west_east'    ), toto, titi & var(0) = titi
-  NCDF_DIMINQ, id, NCDF_DIMID(id, 'south_north'  ), toto, titi & var(1) = titi 
-  NCDF_DIMINQ, id, NCDF_DIMID(id, 'bottom_top'   ), toto, titi & var(2) = titi
-  NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time'         ), toto, titi & var(3) = titi 
-  ;print, var
-
-  ;
-  ; define default data windows
-  ;
-  if (n_elements(offset) le 3) then offset=[0,0,0,0]
-  if (n_elements(stride) le 3) then stride=[1,1,1,1]
-  if (n_elements(count) le 3) then count=intarr(4)
-  w = where(count eq 0) & if (w(0) ne -1) then count[w] = var[w] / stride[w]
-
-  ;
-  ; read variable
-  ;
-  varid=ncdf_varid(id,charvar) ;& print, 'get '+charvar 
-  if (n_elements(notime) eq 0) then timetime = SYSTIME(1) ;& help, /memory
-  ncdf_varget, id, varid, invar, count=count, offset=offset, stride=stride & if (n_elements(notime) eq 0) then print, SYSTIME(1) - timetime, ' s' ;& help, /memory
-  ncdf_close, id
-
-  ;
-  ; OPTIONAL: ANOMALY (value is returned in the keyword)
-  ; 
-  if (n_elements(anomaly) ne 0) then begin
-    anomaly=TOTAL(TOTAL(invar,1),1)/float(var[0])/float(var[1])
-    for i=0,var[0]-1 do for j=0,var[1]-1 do invar(i,j,*,*) = TEMPORARY(invar(i,j,*,*)) - anomaly  ; coute cher ; pas tant
-  endif
-
-  ;
-  ; OPTIONAL: MEAN 
-  ;
-  if (n_elements(coordmean) ne 0) then begin
-  for i=0,nc-1 do begin
-   print, 'mean over dimension '+string(coordmean(i),'(I0)')+' - '+string(var(coordmean(i)),'(I0)')+' elements'
-   invar=total(TEMPORARY(invar),coordmean(i))/var(coordmean(i))
-  endfor
-  endif
-
-  ;
-  ; SAVE
-  ;
-  if (n_elements(nosave) eq 0) then save, invar, filename=name_idl
-
-ENDELSE
-
-if (n_elements(notime) eq 0) then print, SYSTIME(1) - timetime, ' s'
-
-;
-; RETURN FIELD
-;
-return, invar
-
-end
Index: trunk/mesoscale/PLOT/SPEC/MAP/map_uvt.pro
===================================================================
--- trunk/mesoscale/PLOT/SPEC/MAP/map_uvt.pro	(revision 85)
+++ trunk/mesoscale/PLOT/SPEC/MAP/map_uvt.pro	(revision 86)
@@ -24,6 +24,8 @@
 ;
 if (n_elements(field1) ne 0) then getcdf, file=filename, charvar=field1, invar=cfield1
-u = getget(filename, 'Um', count=[0,0,1,1], offset=[0,0,nlevel,ntime])
-v = getget(filename, 'Vm', count=[0,0,1,1], offset=[0,0,nlevel,ntime])
+if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin
+   u = getget(filename, 'Um', count=[0,0,1,1], offset=[0,0,nlevel,ntime])
+   v = getget(filename, 'Vm', count=[0,0,1,1], offset=[0,0,nlevel,ntime])
+endif
 getcdf, file=filename, charvar='XLONG', invar=longi
 getcdf, file=filename, charvar='XLAT', invar=lati
@@ -42,6 +44,8 @@
 	if (no3d ne 'true') then cfield1 = cfield1 (sp:nx-sp-1,sp:ny-sp-1,*,*) else cfield1 = cfield1 (sp:nx-sp-1,sp:ny-sp-1,*)
 endif
-u 		= u		(sp:nx-sp-1,sp:ny-sp-1);,*,*) 
-v       	= v       	(sp:nx-sp-1,sp:ny-sp-1);,*,*)
+if ( (n_elements(overvector_x) ne 0) or (n_elements(field1) eq 0) ) then begin
+  u 		= u		(sp:nx-sp-1,sp:ny-sp-1);,*,*) 
+  v       	= v       	(sp:nx-sp-1,sp:ny-sp-1);,*,*)
+endif
 longi		= longi		(sp:nx-sp-1,sp:ny-sp-1,*)
 lati          	= lati 	        (sp:nx-sp-1,sp:ny-sp-1,*)
@@ -55,6 +59,8 @@
 lon = reform(longi(*,*,ntime))
 lat = reform(lati(*,*,ntime))
-overvector_x = u;reform(u(*,*,nlevel,ntime))
-overvector_y = v;reform(v(*,*,nlevel,ntime))
+if (n_elements(overvector_x) ne 0) then begin
+  overvector_x = u;reform(u(*,*,nlevel,ntime))
+  overvector_y = v;reform(v(*,*,nlevel,ntime))
+endif
 if (n_elements(field1) eq 0) then begin
         print, 'field1: horizontal velocity'
Index: trunk/mesoscale/PLOT/SPEC/getget.pro
===================================================================
--- trunk/mesoscale/PLOT/SPEC/getget.pro	(revision 86)
+++ trunk/mesoscale/PLOT/SPEC/getget.pro	(revision 86)
@@ -0,0 +1,125 @@
+function getget, name, charvar, count=count, offset=offset, stride=stride, coordmean=coordmean, anomaly=anomaly
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; name: char
+; charvar: char
+; count: 4 elements array   ; The dimensions of the sub-image
+; offset: 4 elements array  ; A variable that contains the offset for the sub-image
+; stride: 4 elements array  ; Create a variable to be used as a value for the STRIDE keyword
+; coordmean: 1 to 4 elements array (dimension 1 is 1)
+; anomaly: a non-null element that would be filled with mean
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+
+;print, charvar
+
+
+nosave = 1  ;; comment to get an IDL data file
+notime = 1  ;; comment to get mention of time taken
+if (n_elements(notime) eq 0) then timetime = SYSTIME(1)
+
+;
+; MINIMAL ARGUMENTS
+;
+if (n_elements(name) eq 0) then stop
+if (n_elements(charvar) eq 0) then stop
+
+;
+; IDL DATA DIRECTORY
+;
+name_idl = charvar+name+'.idl'
+  if (n_elements(coordmean) ne 0) then begin
+    nc = n_elements(coordmean)
+    for i=0,nc-1 do begin
+     name_idl = string(coordmean(i),'(I0)')+name_idl   
+    endfor
+  endif
+name_idl = './idl_dat/'+name_idl
+
+;
+; TEST IF FILE EXISTS
+;
+openr,unit,name_idl,/get_lun,error=err
+
+;
+; IF IDL FILE EXISTS, THAT IS RATHER SIMPLE
+;
+IF (err EQ 0) THEN BEGIN
+
+   print, 'read from file'
+   free_lun,unit
+   restore, filename=name_idl
+
+;
+; IF IDL FILE DOES NOT EXIST, WELL, READ NETCDF
+;
+ENDIF ELSE BEGIN
+
+  ;
+  ; check save directory
+  ;
+  if (n_elements(nosave) eq 0) then SPAWN, 'mkdir -p idl_dat'  
+
+  ;
+  ; open file
+  ;
+  id=ncdf_open(name) 
+
+  ;
+  ; get dimensions
+  ;
+  var = intarr(4)
+  NCDF_DIMINQ, id, NCDF_DIMID(id, 'west_east'    ), toto, titi & var(0) = titi
+  NCDF_DIMINQ, id, NCDF_DIMID(id, 'south_north'  ), toto, titi & var(1) = titi 
+  NCDF_DIMINQ, id, NCDF_DIMID(id, 'bottom_top'   ), toto, titi & var(2) = titi
+  NCDF_DIMINQ, id, NCDF_DIMID(id, 'Time'         ), toto, titi & var(3) = titi 
+  ;print, var
+
+  ;
+  ; define default data windows
+  ;
+  if (n_elements(offset) le 3) then offset=[0,0,0,0]
+  if (n_elements(stride) le 3) then stride=[1,1,1,1]
+  if (n_elements(count) le 3) then count=intarr(4)
+  w = where(count eq 0) & if (w(0) ne -1) then count[w] = var[w] / stride[w]
+
+  ;
+  ; read variable
+  ;
+  varid=ncdf_varid(id,charvar) ;& print, 'get '+charvar 
+  if (n_elements(notime) eq 0) then timetime = SYSTIME(1) ;& help, /memory
+  ncdf_varget, id, varid, invar, count=count, offset=offset, stride=stride & if (n_elements(notime) eq 0) then print, SYSTIME(1) - timetime, ' s' ;& help, /memory
+  ncdf_close, id
+
+  ;
+  ; OPTIONAL: ANOMALY (value is returned in the keyword)
+  ; 
+  if (n_elements(anomaly) ne 0) then begin
+    anomaly=TOTAL(TOTAL(invar,1),1)/float(var[0])/float(var[1])
+    for i=0,var[0]-1 do for j=0,var[1]-1 do invar(i,j,*,*) = TEMPORARY(invar(i,j,*,*)) - anomaly  ; coute cher ; pas tant
+  endif
+
+  ;
+  ; OPTIONAL: MEAN 
+  ;
+  if (n_elements(coordmean) ne 0) then begin
+  for i=0,nc-1 do begin
+   print, 'mean over dimension '+string(coordmean(i),'(I0)')+' - '+string(var(coordmean(i)),'(I0)')+' elements'
+   invar=total(TEMPORARY(invar),coordmean(i))/var(coordmean(i))
+  endfor
+  endif
+
+  ;
+  ; SAVE
+  ;
+  if (n_elements(nosave) eq 0) then save, invar, filename=name_idl
+
+ENDELSE
+
+if (n_elements(notime) eq 0) then print, SYSTIME(1) - timetime, ' s'
+
+;
+; RETURN FIELD
+;
+return, invar
+
+end
