Index: /trunk/chantiers/commit_importants.log
===================================================================
--- /trunk/chantiers/commit_importants.log	(revision 107)
+++ /trunk/chantiers/commit_importants.log	(revision 108)
@@ -736,3 +736,70 @@
 Avec la v105, la compilation Titan se fait sans problemes.
 
-
+**-----------
+(commit 106: ajustements dans phytitan)
+**-----------
+- Sorties des opacites (*_hist*.h)
+- g(z) pris en compte (physiq.F, zzlay)
+- correction topo (physiq.F, zzlev(i,1))
+- correction lell.F et lell_light.F (et aussi radtitan.F du coup)
+- remise a jour des sorties
+- suphec.F : RDAY=1.37889e6 et non RSIDAY
+**-----------
+
+*********************
+**** commit_v108 ****
+*********************
+
+Modifs Sebastien Lebonnois dans dyn3d[par] pour sponge layer et dissip horizontale.
+
+** Correction bug
+**---------------
+
+* Modif de dyn3d[par]/calfis[_p].F
+  => ztfi doublement incremente...
+
+** Sponge layer
+**-------------
+
+* Modif de dyn3d[par]/leapfrog[_p].F
+  Utilisation de d*top dans appel a top_bound.F
+  Initialisation des d* (dis, fi, top)
+
+* Modif de dyn3d[par]/top_bound[_p].F
+  Implementation de mode_top_bound
+  Pas d'incrementation
+
+* Modif de dyn3d[par]/conf_gcm.F
+  Parametres de top_bound
+
+* Modif de dyn3d[par]/comconst.h
+  Parametres de top_bound
+
+* Modif de deftank/[titan/venus]/gcm.def
+  Ajout ok_strato
+  Ajout parametres de top_bound
+
+* Ecriture doc
+  Ajout top_bound.tex (et .pdf) dans la documentation.
+
+** Dissipation horizontale
+**------------------------
+
+* Modif de dyn3d[par]/inidissip.F
+  passage a dissip_fac_mid automatique
+  passage a dissip_fac_up sur modele martien si ok_strato
+
+* Modif de dyn3d[par]/conf_gcm.F
+  Parametres dissip_*
+
+* Modif de dyn3d[par]/comconst.h
+  Parametres dissip_*
+
+* Modif de deftank/[titan/venus]/gcm.def
+  Ajout parametres dissip_*   !!! REGLAGES A FAIRE !!!
+
+* Ecriture doc
+  Ajout dissip_horiz.tex (et .pdf) dans la documentation.
+
+Reste la discretisation verticale a voir.
+
Index: /trunk/deftanks/titan/gcm.def
===================================================================
--- /trunk/deftanks/titan/gcm.def	(revision 107)
+++ /trunk/deftanks/titan/gcm.def	(revision 108)
@@ -41,4 +41,26 @@
 ## avec ou sans traceurs                                                 
 iflag_trac=1
+##  Avec ou sans strato // i.e. Couche eponge et second palier pour dissip horiz.
+ok_strato=y
+## Dissipation horizontale
+dissip_fac_mid=2.
+dissip_fac_up=10.
+# deltaz et hdelta en km
+dissip_deltaz=10.
+dissip_hdelta=5.
+# pupstart en Pa
+dissip_pupstart=1.e3
+## Couche eponge 
+#   1: dans les 4 derniers niveaux
+#   2: dans les couches de pression plus faible que 100 fois la pression de la derniere couche
+iflag_top_bound=1
+## Mode Couche eponge 
+#   mode = 0 : pas de sponge
+#   mode = 1 : u et v -> 0
+#   mode = 2 : u et v -> moyenne zonale
+#   mode = 3 : u, v et h -> moyenne zonale
+mode_top_bound=1
+#  Coefficient pour la couche eponge (valeur derniere couche)
+tau_top_bound=4.e-5
 #
 ## longitude en degres du centre du zoom                                 
Index: /trunk/deftanks/venus/gcm.def
===================================================================
--- /trunk/deftanks/venus/gcm.def	(revision 107)
+++ /trunk/deftanks/venus/gcm.def	(revision 108)
@@ -41,4 +41,26 @@
 ## avec ou sans traceurs                                                 
 iflag_trac=0
+##  Avec ou sans strato // i.e. Couche eponge et second palier pour dissip horiz.
+ok_strato=y
+## Dissipation horizontale
+dissip_fac_mid=2.
+dissip_fac_up=2.
+# deltaz et hdelta en km
+dissip_deltaz=10.
+dissip_hdelta=5.
+# pupstart en Pa
+dissip_pupstart=1.e3
+## Couche eponge 
+#   1: dans les 4 derniers niveaux
+#   2: dans les couches de pression plus faible que 100 fois la pression de la derniere couche
+iflag_top_bound=1
+## Mode Couche eponge 
+#   mode = 0 : pas de sponge
+#   mode = 1 : u et v -> 0
+#   mode = 2 : u et v -> moyenne zonale
+#   mode = 3 : u, v et h -> moyenne zonale
+mode_top_bound=1
+#  Coefficient pour la couche eponge (valeur derniere couche)
+tau_top_bound=1.e-5
 #
 ## longitude en degres du centre du zoom                                 
Index: /trunk/documentation/cpdet.tex
===================================================================
--- /trunk/documentation/cpdet.tex	(revision 107)
+++ /trunk/documentation/cpdet.tex	(revision 108)
@@ -112,5 +112,5 @@
 \section{Pratical aspects in the code}
 
-A specific file has been added to the dynamical core, cpdet.F, which includes 
+A specific file has been added to the dynamical core, \textsf{cpdet.F}, which includes 
 all the needed routines to take the $C_p(T)$ possibility into account.
 These routines take advantage of the keyword \textsf{planet\_type} to 
Index: /trunk/documentation/dissip_horiz.tex
===================================================================
--- /trunk/documentation/dissip_horiz.tex	(revision 108)
+++ /trunk/documentation/dissip_horiz.tex	(revision 108)
@@ -0,0 +1,139 @@
+\documentclass[a4paper,10pt]{article}
+%\usepackage{graphicx}
+\usepackage{natbib}  % si appel à bibtex
+%\usepackage[francais]{babel}
+%\usepackage[latin1]{inputenc}  % accents directs (é...), avec babel
+%\usepackage{rotating}
+
+\setlength{\hoffset}{-1.in}
+\setlength{\oddsidemargin}{3.cm}
+\setlength{\textwidth}{15.cm}
+\setlength{\marginparsep}{0.mm}
+\setlength{\marginparwidth}{0.mm}
+
+\setlength{\voffset}{-1.in}
+\setlength{\topmargin}{0.mm}
+\setlength{\headheight}{0.mm}
+\setlength{\headsep}{30.mm}
+\setlength{\textheight}{24.cm}
+\setlength{\footskip}{1.cm}
+
+\setlength{\parindent}{0.mm}
+\setlength{\parskip}{1 em}
+\newcommand{\ten}[1]{$\times 10^{#1}$~} 
+\renewcommand{\baselinestretch}{1.}
+
+\begin{document}
+\pagestyle{plain}
+
+\begin{center}
+{\bf \LARGE 
+Documentation for LMDZ, Planets version
+
+\vspace{1cm}
+\Large
+The horizontal dissipation
+}
+
+\vspace{1cm}
+S\'ebastien Lebonnois
+
+\vspace{1cm}
+Latest version: \today
+\end{center}
+
+\section{Theoretical aspects}
+
+To be written
+
+\section{Pratical aspects in the code}
+
+The horizontal dissipation parameters are chosen in \textsf{gcm.def}.
+
+{\bf Parameters related to the operators}
+
+\begin{itemize}
+\item \textsf{idissip}: timestep for dissipation. 
+Should be equal to \textsf{iperiod}. 
+\item \textsf{lstardis}: boolean that indicates 
+whether to use a star operator (or not). Usually set to {\em True}.
+\item \textsf{nitergdiv}: number of iterations for the {\em gradiv} operator
+\item \textsf{nitergrot}: number of iterations for the {\em nxgradrot} operator
+\item \textsf{niterh}:  number of iterations for the {\em divgrad} operator
+\end{itemize}
+ 
+{\bf Parameters related to timescales}
+
+\begin{itemize}
+\item \textsf{tetagdiv}: time scale (in s) for the {\em gradiv} operator. 
+It corresponds to the attenuation of the smallest wavelengths for u and v 
+perturbations.
+\item \textsf{tetagrot}: time scale (in s) for the {\em nxgradrot} operator. 
+It corresponds to the attenuation of the smallest wavelengths for u and v 
+perturbations.
+\item \textsf{tetatemp}: time scale (in s) for the {\em divgrad} operator. 
+It corresponds to the attenuation of the smallest wavelengths for h 
+perturbations.
+\end{itemize}
+These timescales are the one for the deep atmosphere. 
+However, they are modified by some factors as the pressure decreases. 
+For these factors, two steps are implemented: \textsf{dissip\_fac\_mid} and
+\textsf{dissip\_fac\_up} (also chosen in \textsf{gcm.def}). 
+The first step is always applied, the second is
+applied only when \textsf{ok\_strato} is set to {\em True}.
+Timescales are divided by the factor $f$ computed as detailed below. 
+
+For \textsf{dissip\_fac\_mid} (which is usually equal to 2.), transition is 
+computed with:
+\[
+         x = 1. - \frac{{\rm preff}}{p}
+\]
+\[
+         f_1 = {\rm dissip\_fac\_mid} -
+           \frac{{\rm dissip\_fac\_mid}-1.}{1.+x^2}
+\]
+
+When applied (\textsf{ok\_strato} set to {\em True}), the transition from
+\textsf{dissip\_fac\_mid} to \textsf{dissip\_fac\_up} is done with a 
+$\tanh$ function, using further parameters:
+\begin{itemize}
+\item \textsf{dissip\_deltaz}: altitude range (in km) for the transition.
+\item \textsf{dissip\_hdelta}: scale height (in km) at the altitude of the transition.
+\item \textsf{dissip\_pupstart}: pressure (in Pa) corresponding to the bottom 
+of the transition region.
+\end{itemize}
+The function used is:
+\[
+         f_2 = \left[ 1+
+     \left( \frac{{\rm dissip\_fac\_up}}{{\rm dissip\_fac\_mid}} - 1. \right)
+     \times \left( 1 - \frac{1 + \tanh X}{2} \right) \right]
+\]
+with
+\[
+    X = \frac{6.\times {\rm dissip\_hdelta}}{{\rm dissip\_deltaz}}
+      \log \frac{p}{{\rm dissip\_pupstart} \times 
+     \exp \left( \frac{-{\rm dissip\_deltaz}}{2 {\rm dissip\_hdelta}} \right)}
+%      = \frac{6.\times {\rm dissip\_hdelta}}{{\rm dissip\_deltaz}}
+%      \log \frac{p}{{\rm dissip\_pupstart}} + 3.
+\]
+
+{\bf A last parameter\dots}
+
+There is one last coefficient in \textsf{gcm.def} related to dissipation: 
+\textsf{coefdis}. 
+It is used in \textsf{inigeom.F} to compute \textsf{gamdi\_gdiv, gamdi\_grot, 
+gamdi\_h}. Usually put to 0... Should be explained here...
+
+
+
+
+%\begin{thebibliography}{2}
+%\providecommand{\natexlab}[1]{#1}
+%\expandafter\ifx\csname urlstyle\endcsname\relax
+%  \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
+%  \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
+%  \urlstyle{rm}\Url}\fi
+
+%\end{thebibliography}
+
+\end{document}
Index: /trunk/documentation/top_bound.tex
===================================================================
--- /trunk/documentation/top_bound.tex	(revision 108)
+++ /trunk/documentation/top_bound.tex	(revision 108)
@@ -0,0 +1,88 @@
+\documentclass[a4paper,10pt]{article}
+%\usepackage{graphicx}
+\usepackage{natbib}  % si appel à bibtex
+%\usepackage[francais]{babel}
+%\usepackage[latin1]{inputenc}  % accents directs (é...), avec babel
+%\usepackage{rotating}
+
+\setlength{\hoffset}{-1.in}
+\setlength{\oddsidemargin}{3.cm}
+\setlength{\textwidth}{15.cm}
+\setlength{\marginparsep}{0.mm}
+\setlength{\marginparwidth}{0.mm}
+
+\setlength{\voffset}{-1.in}
+\setlength{\topmargin}{0.mm}
+\setlength{\headheight}{0.mm}
+\setlength{\headsep}{30.mm}
+\setlength{\textheight}{24.cm}
+\setlength{\footskip}{1.cm}
+
+\setlength{\parindent}{0.mm}
+\setlength{\parskip}{1 em}
+\newcommand{\ten}[1]{$\times 10^{#1}$~} 
+\renewcommand{\baselinestretch}{1.}
+
+\begin{document}
+\pagestyle{plain}
+
+\begin{center}
+{\bf \LARGE 
+Documentation for LMDZ, Planets version
+
+\vspace{1cm}
+\Large
+The upper boundary sponge layer
+}
+
+\vspace{1cm}
+S\'ebastien Lebonnois
+
+\vspace{1cm}
+Latest version: \today
+\end{center}
+
+%\section{Theoretical aspects}
+
+\section{Pratical aspects in the code}
+
+The sponge layer is applied at the upper boundary when the \textsf{ok\_strato}
+flag is set to {\em True} in \textsf{gcm.def} 
+(this parameter also controls the application of a second step in the 
+horizontal dissipation).
+
+The tendencies for the upper boundary sponge layer are computed separately in
+the \textsf{top\_bound.F} routine, called from \textsf{leapfrog.F}. 
+These tendencies are \textsf{dutop}, \textsf{dvtop} and \textsf{dhtop}, in
+unit/s. 
+
+Three parameters may be adjusted in the \textsf{gcm.def} file: 
+\begin{itemize}
+\item \textsf{iflag\_top\_bound}: selects the affected layers. 
+  \begin{itemize}
+  \item 1: only the top 4 layers are affected. In this case, the damping rate 
+  is divided by 2 in the second layer, 4 in the third and 8 in the fourth. 
+  \item 2: layers with pressure lower than 100 times the top pressure. 
+  In this case, the damping rate depends linearly on the pressure.
+  \end{itemize}
+\item \textsf{mode\_top\_bound}: selects how the fields are affected.
+  \begin{itemize}
+  \item 0: No sponge layer is applied.
+  \item 1: Zonal and meridional winds are damped to zero.
+  \item 2: Zonal and meridional winds are damped to their zonally averaged value.
+  \item 3: Temperature, zonal and meridional winds are damped to their zonally 
+  averaged value.
+  \end{itemize}
+\item \textsf{tau\_top\_bound}: damping rate (in /s) in the top layer.
+\end{itemize}
+
+%\begin{thebibliography}{2}
+%\providecommand{\natexlab}[1]{#1}
+%\expandafter\ifx\csname urlstyle\endcsname\relax
+%  \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
+%  \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
+%  \urlstyle{rm}\Url}\fi
+
+%\end{thebibliography}
+
+\end{document}
Index: /trunk/libf/dyn3d/calfis.F
===================================================================
--- /trunk/libf/dyn3d/calfis.F	(revision 107)
+++ /trunk/libf/dyn3d/calfis.F	(revision 108)
@@ -587,5 +587,4 @@
 
 ! ADAPTATION GCM POUR CP(T)
-         ztfi=ztfi+zdtfi*dtphys
       call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
 
Index: /trunk/libf/dyn3d/comconst.h
===================================================================
--- /trunk/libf/dyn3d/comconst.h	(revision 107)
+++ /trunk/libf/dyn3d/comconst.h	(revision 108)
@@ -6,9 +6,9 @@
 
       COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
-     &                 iflag_top_bound
+     &                 iflag_top_bound,mode_top_bound
       COMMON/comconstr/dtvr,daysec,                                     &
      & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
-     &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
-     &                   ,tau_top_bound,                                &
+     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
+     & ,dissip_pupstart  ,tau_top_bound,                                &
      & daylen,year_day,molmass, ihf
       COMMON/cpdetvenus/nu_venus,t0_venus
@@ -28,6 +28,7 @@
       REAL g ! (m/s2) gravity
       REAL omeg ! (rad/s) rotation rate of the planet
-      REAL dissip_factz,dissip_deltaz,dissip_zref
-      INTEGER iflag_top_bound
+      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
+      REAL dissip_pupstart
+      INTEGER iflag_top_bound,mode_top_bound
       REAL tau_top_bound
       REAL daylen ! length of solar day, in 'standard' day length
Index: /trunk/libf/dyn3d/conf_gcm.F
===================================================================
--- /trunk/libf/dyn3d/conf_gcm.F	(revision 107)
+++ /trunk/libf/dyn3d/conf_gcm.F	(revision 108)
@@ -303,17 +303,28 @@
 ! Parametres controlant la variation sur la verticale des constantes de
 ! dissipation.
-! Pour le moment actifs uniquement dans la version a 39 niveaux
-! avec ok_strato=y
-
-       dissip_factz=4.
-       dissip_deltaz=10.
-       dissip_zref=30.
-       CALL getin('dissip_factz',dissip_factz )
+! Actifs uniquement avec ok_strato=y
+
+       dissip_fac_mid=2.
+       dissip_fac_up=10.
+       dissip_deltaz=10.! Intervalle (km) pour le changement mid / up
+       dissip_hdelta=5. ! scale height (km) dans la zone de la transition(m)
+       dissip_pupstart=1.e3  ! pression (Pa) au bas la transition mid / up
+       CALL getin('dissip_fac_mid',dissip_fac_mid )
+       CALL getin('dissip_fac_up',dissip_fac_up )
        CALL getin('dissip_deltaz',dissip_deltaz )
-       CALL getin('dissip_zref',dissip_zref )
-
+       CALL getin('dissip_hdelta',dissip_hdelta )
+       CALL getin('dissip_pupstart',dissip_pupstart )
+
+! Parametres controlant la couche eponge
+! Actifs uniquement avec ok_strato=y
+!   mode = 0 : pas de sponge
+!   mode = 1 : u et v -> 0
+!   mode = 2 : u et v -> moyenne zonale
+!   mode = 3 : u, v et h -> moyenne zonale
        iflag_top_bound=1
+       mode_top_bound=3
        tau_top_bound=1.e-5
        CALL getin('iflag_top_bound',iflag_top_bound)
+       CALL getin('mode_top_bound',mode_top_bound)
        CALL getin('tau_top_bound',tau_top_bound)
 
Index: /trunk/libf/dyn3d/inidissip.F
===================================================================
--- /trunk/libf/dyn3d/inidissip.F	(revision 107)
+++ /trunk/libf/dyn3d/inidissip.F	(revision 108)
@@ -31,5 +31,5 @@
       INTEGER l,ij,idum,ii
       REAL tetamin
-      REAL pseudoz
+      REAL Pup
 
       REAL ran1
@@ -177,20 +177,38 @@
 c   --------------------------------------------------
 
-      if (ok_strato .and. llm==39) then
-         do l=1,llm
-            pseudoz=8.*log(preff/presnivs(l))
-            zvert(l)=1+
-     s      (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2.
-     s      *(dissip_factz-1.)
-         enddo 
-      else
-         DO l=1,llm
-            zvert(l)=1.
-         ENDDO
-         fact=2.
-         DO l = 1, llm
-            zz      = 1. - preff/presnivs(l)
-            zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
-         ENDDO
+c First step: going from 1 to dissip_fac_mid (in gcm.def)
+c============
+      DO l=1,llm
+         zz      = 1. - preff/presnivs(l)
+         zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
+      ENDDO
+
+         write(*,*) 'Dissipation : '
+         write(*,*) 'Multiplication de la dissipation en altitude :',
+         write(*,*) '  dissip_fac_mid =', dissip_fac_mid
+
+c Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
+c==========================
+c Utilisation de la fonction tangente hyperbolique pour augmenter
+c arbitrairement la dissipation et donc la stabilite du modele a 
+c partir d'une certaine altitude.
+
+c   Le facteur multiplicatif de basse atmosphere etant deja pris 
+c   en compte, il faut diviser le facteur multiplicatif de haute 
+c   atmosphere par celui-ci.
+
+      if (ok_strato) then
+
+       Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
+       do l=1,llm
+         zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1)
+     &   *(1-(0.5*(1+tanh(-6./dissip_deltaz*
+     &    (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
+       enddo 
+
+         write(*,*) '  dissip_fac_up =', dissip_fac_up
+         write(*,*) 'Transition mid /up:  Pupstart,delta =',
+     &             dissip_pupstart,'Pa', dissip_deltaz , '(km)'
+
       endif
 
Index: /trunk/libf/dyn3d/leapfrog.F
===================================================================
--- /trunk/libf/dyn3d/leapfrog.F	(revision 107)
+++ /trunk/libf/dyn3d/leapfrog.F	(revision 108)
@@ -111,6 +111,4 @@
       REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
       REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
-
-      real :: duspg(ip1jmp1,llm) ! for bilan_dyn
 
 c   variables pour le fichier histoire
@@ -463,14 +461,11 @@
 c      -------------------
          IF (ok_strato) THEN
-           CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
+           CALL top_bound( vcov,ucov,teta,masse,dutop,dvtop,dtetatop)
+c dqtop=0, dptop=0
+           CALL addfi( dtphys, leapf, forward   ,
+     $                  ucov, vcov, teta , q   ,ps ,
+     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
          ENDIF
-c dqtop=0, dptop=0
-       
-c      ajout des tendances physiques:
-c      ------------------------------
-          CALL addfi( dtphys, leapf, forward   ,
-     $                  ucov, vcov, teta , q   ,ps ,
-     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
-c
+
 c  Diagnostique de conservation de l'énergie : difference
          IF (ip_ebil_dyn.ge.1 ) THEN 
@@ -512,10 +507,6 @@
         ! Sponge layer (if any)
         IF (ok_strato) THEN
-          dutop(:,:)=0.
-          dvtop(:,:)=0.
-          dtetatop(:,:)=0.
-          dqtop(:,:,:)=0.
-          dptop(:)=0.
           CALL top_bound(vcov,ucov,teta,masse,dutop,dvtop,dtetatop)
+c dqtop=0, dptop=0
           CALL addfi( dtvr, leapf, forward   ,
      $                  ucov, vcov, teta , q   ,ps ,
@@ -669,5 +660,5 @@
                  CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
      &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
-     &                 du,dudis,duspg,dufi)
+     &                 du,dudis,dutop,dufi)
 #endif
                END IF
@@ -801,5 +792,5 @@
                  CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
      &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
-     &                 du,dudis,duspg,dufi)
+     &                 du,dudis,dutop,dufi)
 #endif
                ENDIF
Index: /trunk/libf/dyn3d/top_bound.F
===================================================================
--- /trunk/libf/dyn3d/top_bound.F	(revision 107)
+++ /trunk/libf/dyn3d/top_bound.F	(revision 108)
@@ -52,4 +52,5 @@
       LOGICAL,SAVE :: first=.true.
 
+      REAL zkm
       INTEGER j,l
 
@@ -68,5 +69,5 @@
              rdamp(llm-3)=tau_top_bound/8.
          else if (iflag_top_bound.eq.2) then
-! couce eponge dans toutes les couches de pression plus faible que
+! couche eponge dans toutes les couches de pression plus faible que
 ! 100 fois la pression de la derniere couche
              rdamp(:)=tau_top_bound
@@ -74,10 +75,26 @@
          endif
          first=.false.
-         print*,'TOP_BOUND rdamp=',rdamp
+         print*,'TOP_BOUND mode',mode_top_bound
+         print*,'Coeffs pour la couche eponge a l equateur'
+         print*,'p (Pa)  z(km)  tau (s)   dt*rdamp'
+         do l=1,llm
+           if (rdamp(l).ne.0.) then
+            zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
+          print*,presnivs(l),zkm,
+     .          1./rdamp(l),
+     .          dt*rdamp(l)
+           endif
+         enddo
       endif
 
       CALL massbar(masse,massebx,masseby)
 
-      do l=1,llm
+c   mode = 0 : pas de sponge
+c   mode = 1 : u et v -> 0
+c   mode = 2 : u et v -> moyenne zonale
+c   mode = 3 : u, v et h -> moyenne zonale
+
+      if (mode_top_bound.ge.2) then
+       do l=1,llm
         do j=1,jjm
           vzon(j,l)=0.
@@ -92,15 +109,7 @@
           vzon(j,l)=vzon(j,l)/zm
         enddo
-      enddo
+       enddo
 
-      do l=1,llm
-        do i=1,iip1
-          do j=1,jjm
-            dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
-          enddo
-        enddo
-      enddo
-
-      do l=1,llm
+       do l=1,llm
         do j=2,jjm
           uzon(j,l)=0.
@@ -112,7 +121,12 @@
           uzon(j,l)=uzon(j,l)/zm
         enddo
-      enddo
+       enddo
+      else
+          vzon(:,:)=0.
+          uzon(:,:)=0.
+      endif
 
-      do l=1,llm
+      if (mode_top_bound.ge.3) then
+       do l=1,llm
         do j=2,jjm
           zm=0.
@@ -124,19 +138,37 @@
           tzon(j,l)=tzon(j,l)/zm
         enddo
-      enddo
+       enddo
+      endif
 
 C   AMORTISSEMENTS LINEAIRES:
 
-      do l=1,llm
+      if (mode_top_bound.ge.1) then
+       do l=1,llm
+        do i=1,iip1
+          do j=1,jjm
+            dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
+          enddo
+        enddo
+       enddo
+
+       do l=1,llm
         do i=1,iip1
           do j=2,jjm
-            du(i,j,l)=du(i,j,l)
-     s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
-            dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
+            du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
           enddo
         enddo
-      enddo
+       enddo
+      endif
+
+      if (mode_top_bound.ge.3) then
+       do l=1,llm
+        do i=1,iip1
+          do j=2,jjm
+            dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
+          enddo
+        enddo
+       enddo
+      endif
       
-
       RETURN
       END
Index: /trunk/libf/dyn3dpar/addfi_p.F
===================================================================
--- /trunk/libf/dyn3dpar/addfi_p.F	(revision 107)
+++ /trunk/libf/dyn3dpar/addfi_p.F	(revision 108)
@@ -28,5 +28,5 @@
 c      pdufi(ip1jmp1,llm)     |
 c      pdvfi(ip1jm,llm)       |   respective
-c      pdhfi(ip1jmp1)         |      tendencies
+c      pdhfi(ip1jmp1)         |      tendencies  (unit/s)
 c      pdtsfi(ip1jmp1)        |
 c
Index: /trunk/libf/dyn3dpar/calfis_p.F
===================================================================
--- /trunk/libf/dyn3dpar/calfis_p.F	(revision 107)
+++ /trunk/libf/dyn3dpar/calfis_p.F	(revision 108)
@@ -952,9 +952,4 @@
       
 ! ADAPTATION GCM POUR CP(T)
-!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-        ztfi(1:klon,l)=ztfi(1:klon,l)+zdtfi(1:klon,l)*dtphys
-      ENDDO
-!$OMP END DO
       call t2tpot_p(ngridmx,llm,ztfi,zteta,zpk)
 
Index: /trunk/libf/dyn3dpar/comconst.h
===================================================================
--- /trunk/libf/dyn3dpar/comconst.h	(revision 107)
+++ /trunk/libf/dyn3dpar/comconst.h	(revision 108)
@@ -6,9 +6,9 @@
 
       COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
-     &                 iflag_top_bound
+     &                 iflag_top_bound,mode_top_bound
       COMMON/comconstr/dtvr,daysec,                                     &
      & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
-     &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
-     &                   ,tau_top_bound,                                &
+     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
+     & ,dissip_pupstart  ,tau_top_bound,                                &
      & daylen,year_day,molmass, ihf
       COMMON/cpdetvenus/nu_venus,t0_venus
@@ -28,6 +28,7 @@
       REAL g ! (m/s2) gravity
       REAL omeg ! (rad/s) rotation rate of the planet
-      REAL dissip_factz,dissip_deltaz,dissip_zref
-      INTEGER iflag_top_bound
+      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
+      REAL dissip_pupstart
+      INTEGER iflag_top_bound,mode_top_bound
       REAL tau_top_bound
       REAL daylen ! length of solar day, in 'standard' day length
Index: /trunk/libf/dyn3dpar/conf_gcm.F
===================================================================
--- /trunk/libf/dyn3dpar/conf_gcm.F	(revision 107)
+++ /trunk/libf/dyn3dpar/conf_gcm.F	(revision 108)
@@ -314,17 +314,28 @@
 ! Parametres controlant la variation sur la verticale des constantes de
 ! dissipation.
-! Pour le moment actifs uniquement dans la version a 39 niveaux
-! avec ok_strato=y
-
-       dissip_factz=4.
-       dissip_deltaz=10.
-       dissip_zref=30.
-       CALL getin('dissip_factz',dissip_factz )
+! Actifs uniquement avec ok_strato=y
+
+       dissip_fac_mid=2.
+       dissip_fac_up=10.
+       dissip_deltaz=10.! Intervalle (km) pour le changement mid / up
+       dissip_hdelta=5. ! scale height (km) dans la zone de la transition(m)
+       dissip_pupstart=1.e3  ! pression (Pa) au bas la transition mid / up
+       CALL getin('dissip_fac_mid',dissip_fac_mid )
+       CALL getin('dissip_fac_up',dissip_fac_up )
        CALL getin('dissip_deltaz',dissip_deltaz )
-       CALL getin('dissip_zref',dissip_zref )
-
+       CALL getin('dissip_hdelta',dissip_hdelta )
+       CALL getin('dissip_pupstart',dissip_pupstart )
+
+! Parametres controlant la couche eponge
+! Actifs uniquement avec ok_strato=y
+!   mode = 0 : pas de sponge
+!   mode = 1 : u et v -> 0
+!   mode = 2 : u et v -> moyenne zonale
+!   mode = 3 : u, v et h -> moyenne zonale
        iflag_top_bound=1
+       mode_top_bound=3
        tau_top_bound=1.e-5
        CALL getin('iflag_top_bound',iflag_top_bound)
+       CALL getin('mode_top_bound',mode_top_bound)
        CALL getin('tau_top_bound',tau_top_bound)
 
Index: /trunk/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- /trunk/libf/dyn3dpar/leapfrog_p.F	(revision 107)
+++ /trunk/libf/dyn3dpar/leapfrog_p.F	(revision 108)
@@ -113,7 +113,9 @@
       REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
 
-       !! Aymeric -- cp(T) comme dans leapfrog.F, SAVE OK ???
-      REAL,SAVE :: duspg(ip1jmp1,llm) ! for bilan_dyn
-
+c   tendances top_bound (sponge layer)
+      REAL,SAVE :: dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
+      REAL,SAVE :: dtetatop(ip1jmp1,llm)
+      REAL,SAVE :: dptop(ip1jmp1)
+      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqtop
 
 c   variables pour le fichier histoire
@@ -242,4 +244,5 @@
             ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
             ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
+            ALLOCATE(dqtop(ip1jmp1,llm,nqtot))
          END IF
 c$OMP END MASTER      
@@ -251,5 +254,20 @@
 
 c$OMP MASTER
+c INITIALISATIONS
+        dudis(:,:)   =0.
+        dvdis(:,:)   =0.
+        dtetadis(:,:)=0.
+        dutop(:,:)   =0.
+        dvtop(:,:)   =0.
+        dtetatop(:,:)=0.
+        dqtop(:,:,:) =0.
+        dptop(:)     =0.
+        dufi(:,:)   =0.
+        dvfi(:,:)   =0.
+        dtetafi(:,:)=0.
+        dqfi(:,:,:) =0.
+        dpfi(:)     =0.
       dq(:,:,:)=0.
+
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
@@ -914,12 +932,18 @@
 c      ajout des tendances physiques:
 c      ------------------------------
-         IF (ok_strato) THEN
-           CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
-         ENDIF
-       
           CALL addfi_p( dtphys, leapf, forward   ,
      $                  ucov, vcov, teta , q   ,ps ,
      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
 
+c      Couche superieure :
+c      -------------------
+         IF (ok_strato) THEN
+           CALL top_bound_p( vcov,ucov,teta,masse,dutop,dvtop,dtetatop)
+           CALL addfi_p( dtphys, leapf, forward   ,
+     $                  ucov, vcov, teta , q   ,ps ,
+     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
+
+         ENDIF
+       
 c$OMP BARRIER
 c$OMP MASTER
@@ -1444,5 +1468,5 @@
                  CALL bilan_dyn_p(dtvr*iperiod,dtvr*day_step*periodav,
      &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
-     &                 du,dudis,duspg,dufi)
+     &                 du,dudis,dutop,dufi)
 c$OMP END MASTER
               ENDIF !ok_dynzon
@@ -1670,5 +1694,5 @@
                  CALL bilan_dyn_p(dtvr*iperiod,dtvr*day_step*periodav,
      &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
-     &                 du,dudis,duspg,dufi)
+     &                 du,dudis,dutop,dufi)
 
 c$OMP END MASTER
Index: /trunk/libf/dyn3dpar/top_bound_p.F
===================================================================
--- /trunk/libf/dyn3dpar/top_bound_p.F	(revision 107)
+++ /trunk/libf/dyn3dpar/top_bound_p.F	(revision 108)
@@ -64,5 +64,5 @@
              rdamp(llm-3)=tau_top_bound/8.
          else if (iflag_top_bound == 2) then
-! couce eponge dans toutes les couches de pression plus faible que
+! couche eponge dans toutes les couches de pression plus faible que
 ! 100 fois la pression de la derniere couche
              rdamp(:)=tau_top_bound
@@ -70,5 +70,15 @@
          endif
          first=.false.
-         print*,'TOP_BOUND rdamp=',rdamp
+         print*,'TOP_BOUND mode',mode_top_bound
+         print*,'Coeffs pour la couche eponge a l equateur'
+         print*,'p (Pa)  z(km)  tau (s)   dt*rdamp'
+         do l=1,llm
+           if (rdamp(l).ne.0.) then
+            zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
+          print*,presnivs(l),zkm,
+     .          1./rdamp(l),
+     .          dt*rdamp(l)
+           endif
+         enddo
 c$OMP END MASTER
 c$OMP BARRIER
@@ -77,4 +87,12 @@
 
       CALL massbar_p(masse,massebx,masseby)
+
+c   mode = 0 : pas de sponge
+c   mode = 1 : u et v -> 0
+c   mode = 2 : u et v -> moyenne zonale
+c   mode = 3 : u, v et h -> moyenne zonale
+
+C POUR V
+
 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
 
@@ -84,5 +102,6 @@
 
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      do l=1,llm
+      if (mode_top_bound.ge.2) then
+       do l=1,llm
         do j=jjb,jje
           zm=0.
@@ -97,16 +116,31 @@
           vzon(j,l)=vzon(j,l)/zm
         enddo
-      enddo
+       enddo
+      else
+       do l=1,llm
+        do j=jjb,jje
+          vzon(j,l)=0.
+        enddo
+       enddo
+      endif
 c$OMP END DO NOWAIT   
 
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      do l=1,llm
+C   AMORTISSEMENTS LINEAIRES:
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+      if (mode_top_bound.ge.1) then
+       do l=1,llm
         do j=jjb,jje
           do i=1,iip1
-            dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
-          enddo
-        enddo
-      enddo
-c$OMP END DO NOWAIT
+            dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
+          enddo
+        enddo
+       enddo
+      endif
+c$OMP END DO NOWAIT
+
+C POUR U ET H
+
+C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
 
       jjb=jj_begin
@@ -116,5 +150,6 @@
 
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      do l=1,llm
+      if (mode_top_bound.ge.2) then
+       do l=1,llm
         do j=jjb,jje
           uzon(j,l)=0.
@@ -126,9 +161,17 @@
           uzon(j,l)=uzon(j,l)/zm
         enddo
-      enddo
+       enddo
+      else
+       do l=1,llm
+        do j=jjb,jje
+          uzon(j,l)=0.
+        enddo
+       enddo
+      endif
 c$OMP END DO NOWAIT
 
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
-      do l=1,llm
+      if (mode_top_bound.ge.3) then
+       do l=1,llm
         do j=jjb,jje
           zm=0.
@@ -140,5 +183,6 @@
           tzon(j,l)=tzon(j,l)/zm
         enddo
-      enddo
+       enddo
+      endif
 c$OMP END DO NOWAIT
 
@@ -146,13 +190,25 @@
 
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      do l=1,llm
+      if (mode_top_bound.ge.1) then
+       do l=1,llm
         do j=jjb,jje
           do i=1,iip1
-            du(i,j,l)=du(i,j,l)
-     s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
-            dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
-          enddo
-       enddo
-      enddo
+            du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
+          enddo
+        enddo
+       enddo
+      endif
+c$OMP END DO NOWAIT
+      
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+      if (mode_top_bound.ge.3) then
+       do l=1,llm
+        do j=jjb,jje
+          do i=1,iip1
+            dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
+          enddo
+        enddo
+       enddo
+      endif
 c$OMP END DO NOWAIT
       
