Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/stratosphere_mask.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/stratosphere_mask.F90	(revision 3024)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/stratosphere_mask.F90	(revision 3356)
@@ -2,5 +2,5 @@
 ! $Id$
 !
-SUBROUTINE stratosphere_mask(missing_val, t_seri, pplay, xlat)
+SUBROUTINE stratosphere_mask(missing_val, pphis, t_seri, pplay, xlat)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -8,23 +8,25 @@
 ! determination of tropopause height and temperature from gridded temperature data
 !
-! reference: Reichler, T., M. Dameris, and R. Sausen (2003):
+! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003)
 ! modified: 6/28/06 tjr
 ! adapted to LMDZ by C. Kleinschmitt (2016-02-15)
-!
-! input:    temp(nlon,nlat,nlev)    3D-temperature field 
-!       ps(nlon,nlat)       2D-surface pressure field
-!       zs(nlon,nlat)       2D-surface height
-!       nlon            grid points in x
-!       nlat            grid points in y
-!       pfull(nlon,nlat,nlev)   full pressure levels in Pa
-!       plimu           upper limit for tropopause pressure
-!       pliml           lower limit for tropopause pressure 
-!       gamma           tropopause criterion, e.g. -0.002 K/m
-!
-! output:   p_tropopause(klon)     tropopause pressure in Pa with missing values
-!           t_tropopause(klon)     tropopause temperature in K with missing values
-!           z_tropopause(klon)     tropopause height in m with missing values
-!           stratomask             stratospheric mask withtout missing values
-!           ifil                  # of undetermined values
+! committed to LMDz by O. Boucher (2016) with a mistake 
+! mistake corrected by O. Boucher (2017-12-11)
+!
+! input:  temp(nlon,nlat,nlev)  3D-temperature field 
+!         ps(nlon,nlat)         2D-surface pressure field
+!         zs(nlon,nlat)         2D-surface height
+!         nlon                  grid points in x
+!         nlat                  grid points in y
+!         pfull(nlon,nlat,nlev) full pressure levels in Pa
+!         plimu                 upper limit for tropopause pressure
+!         pliml                 lower limit for tropopause pressure 
+!         gamma                 tropopause criterion, e.g. -0.002 K/m
+!
+! output: p_tropopause(klon)    tropopause pressure in Pa with missing values
+!         t_tropopause(klon)    tropopause temperature in K with missing values
+!         z_tropopause(klon)    tropopause height in m with missing values
+!         stratomask            stratospheric mask withtout missing values
+!         ifil                  # of undetermined values
 !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -37,5 +39,8 @@
 IMPLICIT NONE
 
+INCLUDE "YOMCST.h"
+
 REAL, INTENT(IN)                       :: missing_val ! missing value, also XIOS
+REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! Geopotentiel de surface
 REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
 REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
@@ -60,5 +65,5 @@
   ENDDO
   psrf=pplay(i,1)
-  zsrf=0.0
+  zsrf=pphis(i)/RG           !--altitude de la surface
   call twmo(missing_val, klev, t, p, psrf, zsrf, plimu, pliml, gamma, ptrp, ttrp, ztrp)
   tp(i)=ptrp
@@ -104,4 +109,6 @@
 subroutine twmo(missing_val, level, t, p, ps, zs, plimu, pliml, gamma, ptrp, ttrp, ztrp)
 
+! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003)
+
 implicit none
 
@@ -128,5 +135,5 @@
 ztrp=missing_val
 
-faktor = -RG/R
+faktor = -RG/RD
 
 do j=level,2,-1
@@ -218,14 +225,14 @@
 TDLNP = TDLNP + TM*DLNP
 
-ZTRP = ZS + TDLNP*R/RG
-
-if (ZTRP .lt. 0) then
-  print*,ZTRP
-  print*,PS
-  print*,P
-  print*,T
-  print*,ZS
-  stop
-endif
+ZTRP = ZS + TDLNP*RD/RG
+
+!!if (ZTRP .lt. 0) then
+!!  print*,'ZTRP=',ZTRP
+!!  print*,'PS=',PS
+!!  print*,'P=',P
+!!  print*,'T=',T
+!!  print*,'ZS=',ZS
+!!  stop
+!!endif
 
 return
