Index: LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/physiq_mod.F90	(revision 2979)
+++ LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/physiq_mod.F90	(revision 2980)
@@ -195,6 +195,5 @@
        beta_prec,  &
        rneb,  &
-       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
-       pr_tropopause
+       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
        !
     USE phys_state_var_mod ! Variables sauvegardees de la physique
@@ -2017,9 +2016,8 @@
           IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
           IF(adjust_tropopause) THEN
-             CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot, &
-                  pr_tropopause)
              CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
                       ro3i,  press_edg_climoz,  paprs,   wo, time_climoz,    &
-                           latitude_deg,  press_cen_climoz,  pr_tropopause)
+                      longitude_deg,   latitude_deg,    press_cen_climoz,    &
+                      dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
           ELSE
              CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
Index: LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/regr_pr_time_av_m.F90
===================================================================
--- LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/regr_pr_time_av_m.F90	(revision 2979)
+++ LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/regr_pr_time_av_m.F90	(revision 2980)
@@ -2,7 +2,7 @@
 MODULE regr_pr_time_av_m
 
-  USE print_control_mod, ONLY: lunout
   USE write_field_phy
   USE mod_phys_lmdz_transfert_para, ONLY: bcast
+  USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank
   IMPLICIT NONE
 
@@ -14,12 +14,10 @@
 !  * in all the threads of all the processes, regrid the fields in pressure
 !    to the LMDZ vertical grid.
-!  * input files fields are stretched in the viscinity (+/- 5 kms) of the mean
-!    tropopause (geometrical mean of LMDZ and input fields tropopauses) to force
-!    the resulting field tropopause height to the one of LMDZ. To switch this
-!    feature on, the following arguments must be present:
+!  * the forcing fields are stretched if the following arguments are present:
 !     - "lat_in":   input file latitudes.
 !     - "pcen_in":  input file cells center pressure levels.
 !     - "ptrop_ou": target grid (LMDZ) tropopause pressure.
-!    Note that the ozone quantity conservation is not ensured.
+!   so that the tropopause is shifted to the position of the one of LMDZ.
+!  Note that the ozone quantity conservation is not ensured.
 !-------------------------------------------------------------------------------
 ! Initial routine: regr_pr_av_m module (L. Guez).
@@ -27,6 +25,6 @@
 !    - time interpolation
 !    - 3D compliant
-!    - vertical stretching of the field to allow tropopause and ground pressure
-!      matching between input field and current lmdz field.
+!    - vertical stretching of the field to allow tropopause matching between
+!    input field and current lmdz field.
 !-------------------------------------------------------------------------------
 ! Remarks:
@@ -43,19 +41,32 @@
 !  * The input file pressure at tropopause can be (in decreasing priority):
 !    1) read from the file if "tropopause_air_pressure" field is available.
-!    2) computed from the input file ozone field using:
-!       - o3 concentration at tropopause if "tro3_at_tropopause" is available.
-!       - a default value (100ppbv) if not.
+!    2) computed using "tro3" and "tro3_at_tropopause' (if available).
+!    3) computed using "tro3" and a fixed threshold otherwise, determined using
+!    an empirical three parameters law:
+!         o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3)
+!       => co1 and co2 are in ppbV, and co3 in degrees.
+!       => co3 allow a smooth transition between north and south hemispheres.
 !  * If available, the field "ps" (input file pressure at surface) is used.
 !    If not, the current LMDZ ground pressure is taken instead.
-!  * The O3 threshold for tropopause is defined using 3 coefficients:
-!    o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/c)
-!     => co1 and co2 are in ppmV, and co3 in degrees.
-!     => co3 allow a smooth transition between north and south hemispheres.
-!-------------------------------------------------------------------------------
-! * Fields with suffix "m"/"p" are at the closest records earlier/later than the
-!   middle of the julian day "julien", on the global "dynamics" horizontal grid.
-! * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
-!   interval "pint_in(k:k+1)]" and variable "nam(l)".
-! * In the 2D file case, the values are the same for all longitudes.
+!  * Fields with suffix "m"/"p" are at the closest records earlier/later than
+!  the mid-julian day "julien", on the global "dynamics" horizontal grid.
+!  * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
+!    interval "pint_in(k:k+1)]" and variable "nam(l)".
+!  * In the 2D file case, the values are the same for all longitudes.
+!  * The tropopause correction works like this: the input fields (file) are
+!  interpolated on output (LMDZ) pressure field, which is streched using a power
+!  law in a limited zone made of 3 layers:
+!    1) between the two tropopauses (file and LMDZ)
+!    2) in an upper and a lower transitions layers
+!  The stretching function has the following form:
+!    Sigma_str = Sigma^(1+alpha), with alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
+!  This value shifts the file tropopause to the height of the one of LMDZ.
+!  The stretching is fully applied in the central zone only, and only partially
+!  in the transitions zones, thick enough to guarantee a growing stretched
+!  pressure field. The ponderation function for alpha to modulate the stretching
+!  is constant equal to 1 in the central layer, and quasi-linear (from 1 to 0)
+!  in the transition layers (from 1 to 0 ; in fact: sections of 1/log function),
+!  making this localisation function quasi-trapezoidal.
+!
 ! * The following fields are on the global "dynamics" grid, as read from files:
   REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
@@ -64,4 +75,6 @@
   REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
   REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
+  INTEGER, SAVE :: ntim_in                         !--- Records nb in input file
+  INTEGER, SAVE :: itrp0                           !--- idx above chem tropop.
   INTEGER, SAVE :: irec                            !--- Current time index
 !      * for daily   input files: current julian day number
@@ -73,7 +86,9 @@
   LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
 !$OMP THREADPRIVATE(lfirst)
-  REAL, PARAMETER :: co3(3)=[91.,28.,20.]          !--- Coeffs for o3 threshold
-  REAL, PARAMETER :: prtrop=5.E+3  !--- Value lower than the tropopause pressure
-  REAL, PARAMETER :: delta =5.     !--- Dist. around tropopause for strain (kms)
+  REAL,    PARAMETER :: pTropUp=5.E+3 !--- Value < tropopause pressure (Pa)
+  REAL,    PARAMETER :: gamm = 0.4    !--- Relative thickness of transitions
+  REAL,    PARAMETER :: rho  = 1.3    !--- Max tropopauses sigma ratio
+  REAL,    PARAMETER :: o3t0 = 1.E-7  !--- Nominal O3 vmr at tropopause
+  LOGICAL, PARAMETER :: lo3tp=.FALSE. !--- Use paramztrized O3 vmr at tropopause
 
 CONTAINS
@@ -81,10 +96,11 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3, time_in,    &
-                                     lat_in, pcen_in, ptrop_ou)
+SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3,             &
+                             time_in, lon_in, lat_in, pcen_in, ptrop_ou)
 !
 !-------------------------------------------------------------------------------
   USE dimphy,         ONLY: klon
-  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err
+  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, &
+                            NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION
   USE netcdf,         ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR
   USE assert_m,       ONLY: assert
@@ -105,11 +121,12 @@
   REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
   REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces file Pres, Pa, ascending
-  REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ Pres, Pa (klon,llm+1)
+  REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ pressure, Pa (g,nbp_lev+1)
   REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
   REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
                                               !    since Jan 1 of current year
+  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- Input/output longitudes vector
   REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
   REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure
-  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pressure (klon)
+  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pres (klon)
 !-------------------------------------------------------------------------------
 ! Local variables:
@@ -117,14 +134,18 @@
   include "YOMCST.h"
   CHARACTER(LEN=80)  :: sub
-  CHARACTER(LEN=320) :: msg
-  INTEGER :: vID, ncerr, n_var, nlev_in,ntim_in, ndim, i, ibot, itop, itrp,itrp0
+  CHARACTER(LEN=320) :: str
+  INTEGER :: vID, ncerr, n_var, ibot, ibo0, nn, itrp
+  INTEGER :: i, nlev_in, n_dim, itop, ito0, i0
   LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
   REAL    :: y_frac                           !--- Elapsed year fraction
-  REAL    :: alpha, beta, al                  !--- For strectching/interpolation
-  REAL    :: SigT_in, SigT_ou, SigT_mn        !--- Tropopause: in/out/mean
-  REAL    :: SigA_bot, SigA_top               !--- Strained domain bounds
-  REAL    :: Sig_in (SIZE(pint_in))           !--- Input field sigma levels
-  REAL    :: phi    (SIZE(pint_in))           !--- Stretching exponent anomaly
-  REAl    :: Pint_st(SIZE(pint_in))           !--- Stretched pressure levels
+  REAL    :: alpha, beta, al                  !--- For stretching/interpolation
+  REAL    :: SigT_in, SigT_ou                 !--- Input and output tropopauses
+  REAL    :: Sig_bot, Sig_top                 !--- Fully strained layer  bounds
+  REAL    :: Sig_bo0, Sig_to0                 !--- Total strained layers bounds
+  REAL    :: Sig_in(SIZE(pint_in))            !--- Input field sigma levels
+  REAL    :: Sig_ou (nbp_lev+1)               !--- Output LMDZ sigma levels
+  REAL    :: phi    (nbp_lev+1)               !--- Stretching exponent anomaly
+  REAL    :: pstr_ou(nbp_lev+1)               !--- Stretched pressure levels
+  REAL    :: pintou (nbp_lev+1)               !--- pint_ou in reversed order
   REAL    :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam))
   REAL    :: v2(klon,           SIZE(pint_in)-1,SIZE(nam))
@@ -134,16 +155,19 @@
 !     "v2(i,    k, l)" is at longitude-latitude  "xlon(i)-xlat(i)".
 ! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)"
-  REAL, DIMENSION(nbp_lon, nbp_lat) ::  ps1, pt1, ot1
-  REAL, DIMENSION(klon)             ::  ps2, pt2, ot2
+  REAL, DIMENSION(nbp_lon, nbp_lat)   :: ps1, pt1, ot1
+  REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
+  LOGICAL :: ll
 !-------------------------------------------------------------------------------
   sub="regr_pr_time_av"
-  nlev_in=SIZE(pint_in)-1; ntim_in=SIZE(time_in)
-  CALL assert(SIZE(v3,1)==klon,                 TRIM(sub)//" v3 klon")
-  CALL assert(SIZE(v3,2)==nbp_lev,              TRIM(sub)//" v3 nbp_lev")
-  n_var = assert_eq(SIZE(nam), SIZE(v3,3),      TRIM(sub)//" v3 n_var")
-  CALL assert(SHAPE(pint_ou)==[klon, nbp_lev+1],TRIM(sub)//" pint_ou")
+  nlev_in=SIZE(pint_in)-1
+  CALL assert(SIZE(v3,1)==klon,          TRIM(sub)//" v3 klon")
+  CALL assert(SIZE(v3,2)==nbp_lev,       TRIM(sub)//" v3 nbp_lev")
+  n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
+  CALL assert(SIZE(pint_ou,1)==klon     ,TRIM(sub)//" pint_ou klon")
+  CALL assert(SIZE(pint_ou,2)==nbp_lev+1,TRIM(sub)//" pint_ou nbp_lev+1")
+  IF(PRESENT(lon_in))   CALL assert(SIZE(lon_in  )==klon,TRIM(sub)//" lon_in klon")
   IF(PRESENT(lat_in))   CALL assert(SIZE(lat_in  )==klon,TRIM(sub)//" lat_in klon")
   IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon")
-  IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in)==nlev_in,TRIM(sub)//" pcen_in")
+  IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in )==nlev_in,TRIM(sub)//" pcen_in")
   lAdjTro=PRESENT(ptrop_ou)
   IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) &
@@ -158,5 +182,7 @@
       lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
       lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
-      linterp=PRESENT(time_in); IF(linterp) linterp=ntim_in==14
+      CALL NF95_INQ_DIMID(fID,"time",vID)
+      CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
+      linterp=PRESENT(time_in).AND.ntim_in==14
       IF(linterp) THEN
         ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
@@ -167,16 +193,14 @@
       END IF
       !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
-      IF(lAdjTro) itrp0=locate(pcen_in,prtrop)
-      IF(lPrSurf) WRITE(lunout,*)TRIM(sub)//': Using GROUND PRESSURE from input O3 forcing file.'
-      IF(linterp) WRITE(lunout,*)TRIM(sub)//': Monthly O3 files => ONLINE TIME INTERPOLATION.'
-      IF(lAdjTro) THEN
-        WRITE(lunout,*)TRIM(sub)//': o3 forcing file tropopause location uses:'
-        IF(lPrTrop) THEN
-          WRITE(lunout,*)'    PRESSURE AT TROPOPAUSE from file.'
-        ELSE IF(lO3Trop) THEN
-          WRITE(lunout,*)'    O3 CONCENTRATION AT TROPOPAUSE from file.'
-        ELSE
-          WRITE(lunout,*)'    PARAMETRIZED O3 concentration at tropopause.'
-        END IF
+      IF(lAdjTro) itrp0=locate(pcen_in,pTropUp)
+      CALL msg(lPrSurf,'Using GROUND PRESSURE from input O3 forcing file.',sub)
+      CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.',sub)
+      CALL msg(lAdjTro,'o3 forcing file tropopause location uses:',sub)
+      IF(lPrTrop) THEN
+        CALL msg(lAdjTro,'    PRESSURE AT TROPOPAUSE from file.')
+      ELSE IF(lO3Trop) THEN
+        CALL msg(lAdjTro,'    O3 CONCENTRATION AT TROPOPAUSE from file.')
+      ELSE
+        CALL msg(lAdjTro,'    PARAMETRIZED O3 concentration at tropopause.')
       END IF
     END IF
@@ -185,9 +209,9 @@
     CALL update_fields()
 
-
     !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
     IF(linterp) THEN
-      WRITE(lunout,'(3(a,f7.3))')TRIM(sub)//': Interpolating O3 at julian day '&
-      &,julien,' from fields at times ',time_in(irec),' and ', time_in(irec+1)
+      WRITE(str,'(3(a,f7.3))')'Interpolating O3 at julian day ',julien,' from '&
+      &//'fields at times ',time_in(irec),' and ', time_in(irec+1)
+      CALL msg(.TRUE.,str,sub)
       al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
       v1=al*v1m+(1.-al)*v1p
@@ -199,11 +223,12 @@
   !$OMP END MASTER
   IF(lfirst) THEN
-    lfirst=.FALSE. ;     CALL bcast(lfirst) ; CALL bcast(lPrSurf)
-    CALL bcast(lPrTrop); CALL bcast(lO3Trop); CALL bcast(linterp)
-    IF(lAdjTro) CALL bcast(itrp0)
+    lfirst=.FALSE.;      CALL bcast(lfirst)
+    IF(lAdjTro)          CALL bcast(itrp0)
+    CALL bcast(lPrTrop); CALL bcast(lPrSurf)
+    CALL bcast(lO3Trop); CALL bcast(linterp)
   END IF
   CALL scatter2d(v1,v2)
   !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
-  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=Pint_ou(:,1); END IF
+  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF
   IF(lPrTrop) CALL scatter2d(pt1,pt2)
   IF(lO3Trop) CALL scatter2d(ot1,ot2)
@@ -212,35 +237,75 @@
   IF(.NOT.lAdjTro) THEN
     DO i=1,klon
-      CALL regr_conserv(1,v2(i,:,:) , Pint_in , Pint_ou(i,nbp_lev+1:1:-1),    &
-                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in))
+      pintou = pint_ou(i,nbp_lev+1:1:-1)
+      CALL regr_conserv(1,v2(i,:,:), pint_in(:), pintou(:),                    &
+                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:), pint_in(:)))
     END DO
   ELSE
     y_frac=(REAL(days_elapsed)+jH_cur)/year_len
+
+    !--- OUTPUT SIGMA LEVELS
     DO i=1,klon
-      SigT_in = get_SigTrop(i,itrp)        !--- input  (file) tropopause
-      SigT_ou = ptrop_ou(i)/pint_ou(i,1)   !--- output (lmdz) tropopause
-      SigT_mn = SQRT(SigT_in*SigT_ou)      !--- mean tropopause>strained domain
-
-      !--- DEFINE THE VERTICAL LAYER IN WHICH THE PROFILE IS STRECHED
-      beta=EXP(delta/scaleheight); Sig_in(:) = [pint_in(1:nlev_in)/ps2(i),1.]
-      SigA_bot = SigT_mn*beta ; ibot=locate(Sig_in(:),SigA_bot)
-      SigA_top = SigT_mn/beta ; itop=locate(Sig_in(:),SigA_top)+1
-
-      !--- HAT FUNCTION phi (/=0 in [SigA_bot,SigA_top] only)
+
+      !--- LOCAL INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS AT INTERFACES
+      pintou=pint_ou(i,nbp_lev+1:1:-1)            !--- increasing values
+      Sig_in(:) = [pint_in(1:nlev_in+1)/ps2(i)]   !--- increasing values
+      Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values
+
+      !--- INPUT (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE
+      SigT_in = get_SigTrop(i,itrp)
+      SigT_ou = ptrop_ou(i)/ps2(i)
+
+      !--- AVOID THE FILAMENTS WHICH WOULD NEED A VERY THICK STRETCHED REGION
+      IF(SigT_ou>SigT_in*rho) SigT_ou = SigT_in*rho
+      IF(SigT_ou<SigT_in/rho) SigT_ou = SigT_in/rho
+      ptropou(i)=SigT_ou*ps2(i)
+
+      !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW
+      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
+
+      !--- FULLY STRETCHED LAYER BOUNDS (FILE AND MODEL TROPOPAUSES)
+      Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot)
+      Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top)
+
+      !--- PARTIALLY STRETCHED LAYER BOUNDS, ENSURING >0 DERIVATIVE
+      beta = LOG(Sig_top)/LOG(Sig_bot)
+      Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
+      Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top **  beta
+
+      !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS
+      !--- gamma<log(Sig_bo0/|alpha|) to keep Sig_bo0<1
+      Sig_bo0 = MIN(Sig_bo0*EXP( gamm*ABS(alpha)), 0.95+(1.-0.95)*Sig_bo0)
+      Sig_to0 =     Sig_to0*EXP(-gamm*ABS(alpha))
+      ibo0 = locate(Sig_ou(:),Sig_bo0)
+      ito0 = locate(Sig_ou(:),Sig_to0)
+
+      !--- FUNCTION FOR STRETCHING LOCALISATION
+      !    0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1
       phi(:)=0.
-      phi(itop:itrp)=(Sig_in(itop:itrp)-SigA_top)/(SigT_in-SigA_top)
-      phi(itrp:ibot)=(SigA_bot-Sig_in(itrp:ibot))/(SigA_bot-SigT_in)
-
-      !--- STRAINED INPUT logP PROFILE ; alpha: MAX. STRETCH. EXPONENT INCREMENT
-      alpha = LOG(SigT_ou/SigT_in)/LOG(SigT_in)
-      Pint_st(:) = pint_ou(i,1) * Sig_in(:)**(1+alpha*phi(:))
-
-      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL LEVELS
-      CALL regr_conserv(1,v2(i,:,:) , Pint_st,  Pint_ou(i,nbp_lev+1:1:-1),     &
-                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),Pint_st))
+      phi(itop+1:ibot) =  1.
+      phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))&
+                            *LOG(Sig_top)/LOG(Sig_top/Sig_to0)
+      phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))&
+                            *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)
+
+      !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
+      pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
+
+      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
+      CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:),                  &
+                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:)))
+
+      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
+      i0=nbp_lev-locate(pintou(:),ptropou(i))+1
+      ll=check_ozone(v3(i, 1:i0     ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &
+                     5.E-9,3.0E-7)
+!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
+      ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', &
+                     5.E-8,1.5E-5)
+!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
+
     END DO
   END IF
 
-
 CONTAINS
 
@@ -252,5 +317,5 @@
 !-------------------------------------------------------------------------------
   IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
-    WRITE(lunout,*)TRIM(sub)//': Updating Ozone forcing field: read from file.'
+    CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
     irec=MIN(INT(julien)+1,ntim_in)     !--- irec is just the julian day number
     !--- MIN -> Security in the unlikely case of roundup errors.
@@ -263,5 +328,5 @@
   ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
     IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
-    WRITE(lunout,*)TRIM(sub)//': Refreshing adjacent Ozone forcing fields.'
+    CALL msg(.TRUE.,sub,'Refreshing adjacent Ozone forcing fields.')
     IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
       irec=locate(time_in,julien)       !--- Need to locate surrounding times
@@ -309,11 +374,11 @@
 !-------------------------------------------------------------------------------
   CALL NF95_INQ_VARID(fID, TRIM(var), vID)
-  CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=ndim)
-  IF(ndim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
-  IF(ndim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
+  CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
+  IF(n_dim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
+  IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
   CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
 
   !--- Flip latitudes: ascending in input file, descending in "rlatu".
-  IF(ndim==3) THEN
+  IF(n_dim==3) THEN
     v(1,:) = v(1,nbp_lat:1:-1)
     v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
@@ -341,12 +406,12 @@
   DO i=1,SIZE(nam)
     CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
-    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=ndim)
-    IF(ndim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
-    IF(ndim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
+    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
+    IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
+    IF(n_dim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
     CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
   END DO
 
   !--- Flip latitudes: ascending in input file, descending in "rlatu".
-  IF(ndim==3) THEN
+  IF(n_dim==3) THEN
     v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
     v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
@@ -360,4 +425,5 @@
 
 
+
 !-------------------------------------------------------------------------------
 !
@@ -370,30 +436,180 @@
   REAL                 :: get_Sigtrop
 !-------------------------------------------------------------------------------
+  !--- Pressure at tropopause is read in the forcing file
+  IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
+    get_SigTrop=pt2(ih)/ps2(ih); RETURN
+  END IF
+  !--- Chemical tropopause definition is used using a particular threshold
+  IF(lO3Trop) THEN                             !--- o3trop KNOWN FROM FILE
+    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,ot2(ih))
+  ELSE IF(lo3tp) THEN                          !--- o3trop PARAMETRIZATION
+    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in)
+  ELSE                                         !--- o3trop CONSTANT
+    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,o3t0)
+  END IF
+  get_SigTrop=get_SigTrop/ps2(ih)
+
+END FUNCTION get_SigTrop
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Determine the tropopause using chemical criterium, ie the pressure
+!          at which the ozone concentration reaches a certain level.
+!-------------------------------------------------------------------------------
+! Remarks:
+! * Input field is upside down (increasing pressure // increasing vertical idx)
+!   The sweep is done from top to ground, starting at the 50hPa layer (idx it0),
+!   where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached:
+!   the O3 profile in this zone is decreasing with pressure.
+! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or
+!   determined using an empirical law roughly derived from ... & al.
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL ::                   chem_tropopause    !--- Pressure at tropopause
+  INTEGER,        INTENT(IN)  :: ih            !--- Horizontal index
+  INTEGER,        INTENT(OUT) :: it            !--- Index of tropopause layer
+  INTEGER,        INTENT(IN)  :: it0           !--- Idx: higher than tropopause
+  REAL,           INTENT(IN)  :: pint(:)       !--- Cells-interf Pr, increasing
+  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)
+  REAL, OPTIONAL, INTENT(IN)  :: pcen(:)       !--- Cells-center Pr, increasing
+  REAL, OPTIONAL, INTENT(IN)  :: o3trop        !--- Ozone at tropopause
+!-------------------------------------------------------------------------------
 ! Local variables:
-  INTEGER :: ii                                !--- Idx of layer containing o3t
   REAL    :: o3t                               !--- Ozone concent. at tropopause
-  REAL    :: prt                               !--- Air pressure   at tropopause
   REAL    :: al                                !--- Interpolation coefficient
-  REAL    :: coef                              !--- Coef: North/South transition
-!-------------------------------------------------------------------------------
-  !--- DETERMINE PRESSURE AT TROPOPAUSE AND DIVIDE IT BY GROUND PRESSURE
-  IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
-    get_SigTrop=pt2(ih)/ps2(ih)
-    it=locate(pint_in(:),pt2(ih))
-  ELSE                                         !--- O3 THRESHOLD
-    coef = TANH(lat_in(i)/co3(3))              !--- Latitude dependant ponderat.
-    coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Time-dependant ponderation
+  REAL    :: coef                              !--- Coeff of latitude modulation
+  REAL, PARAMETER :: co3(3)=[91.,28.,20.]      !--- Coeff for o3 at tropopause
+!-------------------------------------------------------------------------------
+  !--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN
+  IF(PRESENT(o3trop)) THEN                     !=== THRESHOLD FROM ARGUMENTS
+    o3t=o3trop
+  ELSE                                         !=== EMPIRICAL LAW
+    coef = TANH(lat_in(ih)/co3(3))             !--- Latitude  modulation
+    coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Seasonal  modulation
     o3t  = 1.E-9 * (co3(1)+co3(2)*coef)        !--- Value from parametrization
-    IF(lO3Trop) o3t=ot2(ih)                    !--- Value from file
-    !--- Starts from 50hPa and go down.
-    it=itrp0; DO WHILE(v2(ih,it+1,1)>=o3t); it=it+1; END DO
-    al=(v2(ih,it,1)-o3t)/(v2(ih,it,1)-v2(ih,it+1,1))
-    get_SigTrop = ( pcen_in(it)**(1.-al) * pcen_in(it+1)**al )/ps2(ih)
-  END IF
-
-END FUNCTION get_SigTrop
-!
-!-------------------------------------------------------------------------------
-
+  END IF
+
+  !--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)
+  it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
+  al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
+  IF(PRESENT(pcen)) THEN
+    chem_tropopause =       pcen(it)**(1.-al) * pcen(it+1)**al
+  ELSE
+    chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) )
+  END IF
+  it = locate(pint(:), chem_tropopause)        !--- pint(it)<ptrop<pint(it+1)
+
+END FUNCTION chem_tropopause
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  LOGICAL                      :: check_ozone      !--- .T. => some wrong values
+  REAL,             INTENT(IN) :: o3col(:), lon, lat
+  INTEGER,          INTENT(IN) :: ilev0
+  CHARACTER(LEN=*), INTENT(IN) :: layer
+  REAL, OPTIONAL,   INTENT(IN) :: vmin
+  REAL, OPTIONAL,   INTENT(IN) :: vmax
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: k
+  LOGICAL :: lmin, lmax
+  REAL    :: fac
+  CHARACTER(LEN=6) :: unt
+!-------------------------------------------------------------------------------
+  !--- NOTHING TO DO
+  lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
+  lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
+  check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN
+
+  !--- SOME TOO LOW VALUES FOUND
+  IF(lmin) THEN
+    CALL unitc(vmin,unt,fac)
+    DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE
+      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
+        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', &
+        fac*vmin,unt//' in '//TRIM(layer)
+    END DO
+  END IF
+
+  !--- SOME TOO HIGH VALUES FOUND
+  IF(lmax) THEN
+    CALL unitc(vmax,unt,fac)
+    DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE
+      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
+        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', &
+        fac*vmax,unt//' in '//TRIM(layer)
+    END DO
+  END IF
+
+END FUNCTION check_ozone
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE unitc(val,unt,fac)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL,             INTENT(IN)  :: val
+  CHARACTER(LEN=6), INTENT(OUT) :: unt
+  REAL,             INTENT(OUT) :: fac
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: ndgt
+!-------------------------------------------------------------------------------
+  ndgt=3*FLOOR(LOG10(val)/3.)
+  SELECT CASE(ndgt)
+    CASE(-6);     unt=' ppmV '; fac=1.E6
+    CASE(-9);     unt=' ppbV '; fac=1.E9
+    CASE DEFAULT; unt=' vmr  '; fac=1.0
+  END SELECT
+
+END SUBROUTINE unitc
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE msg(ll,str,sub)
+!
+!-------------------------------------------------------------------------------
+  USE print_control_mod, ONLY: lunout
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  LOGICAL,                    INTENT(IN) :: ll
+  CHARACTER(LEN=*),           INTENT(IN) :: str
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
+!-------------------------------------------------------------------------------
+  IF(.NOT.ll) RETURN
+  IF(PRESENT(sub)) THEN
+    WRITE(lunout,*)TRIM(sub)//': '//TRIM(str)
+  ELSE
+    WRITE(lunout,*)TRIM(str)
+  END IF
+
+END SUBROUTINE msg
+!
+!-------------------------------------------------------------------------------
 
 END SUBROUTINE regr_pr_time_av
Index: LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/tropopause_m.F90
===================================================================
--- LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/tropopause_m.F90	(revision 2979)
+++ LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/tropopause_m.F90	(revision 2980)
@@ -2,4 +2,5 @@
 
 !  USE phys_local_var_mod, ONLY: ptrop => pr_tropopause
+  IMPLICIT NONE
   PRIVATE
   PUBLIC :: dyn_tropopause
@@ -9,32 +10,35 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE dyn_tropopause(t,ts,paprs,pplay,presnivs,rot,ptrop,tpot0,pvor0,pmin0,pmax0)
+FUNCTION dyn_tropopause(t, ts, paprs, pplay, rot, thet0, pvor0)
 !
 !-------------------------------------------------------------------------------
-  USE assert_m,      ONLY: assert
-  USE assert_eq_m,   ONLY: assert_eq
-  USE comvert_mod,   ONLY: preff
-  USE dimphy,        ONLY: klon, klev
-  USE geometry_mod,  ONLY: latitude_deg, longitude_deg
-  USE interpolation, ONLY: locate
-  IMPLICIT NONE
+  USE assert_m,     ONLY: assert
+  USE assert_eq_m,  ONLY: assert_eq
+  USE dimphy,       ONLY: klon, klev
+  USE geometry_mod, ONLY: latitude_deg, longitude_deg
+  USE vertical_layers_mod, ONLY: aps, bps, preff
+
 !-------------------------------------------------------------------------------
 ! Arguments:
+  REAL ::     dyn_tropopause(klon) !--- Pressure at tropopause
   REAL, INTENT(IN)  ::      t(:,:) !--- Cells-centers temperature
   REAL, INTENT(IN)  ::     ts(:)   !--- Surface       temperature
   REAL, INTENT(IN)  ::  paprs(:,:) !--- Cells-edges   pressure
   REAL, INTENT(IN)  ::  pplay(:,:) !--- Cells-centers pressure
-  REAL, INTENT(IN)  :: presnivs(:) !--- Cells-centers nominal pressure
   REAL, INTENT(IN)  ::    rot(:,:) !--- Cells-centers relative vorticity
-  REAL, INTENT(OUT) :: ptrop(klon) !--- Tropopause pressure
-  REAL, INTENT(IN), OPTIONAL :: tpot0, pvor0, pmin0, pmax0
+  REAL, INTENT(IN), OPTIONAL :: thet0, pvor0
 !-------------------------------------------------------------------------------
 ! Local variables:
   include "YOMCST.h"
   CHARACTER(LEN=80)  :: sub
-  INTEGER :: it, i, nsrf, k, nh, kmin, kmax
-  REAL    :: p1, p2, tpo0, pvo0, pmn0, pmx0, al
-  REAL, DIMENSION(klon,klev)   :: t_edg, t_pot
-  REAL, DIMENSION(klon,klev-1) :: pvort
+  INTEGER :: i, k, kb, kt, kp, ib, ie, nw
+  REAL    :: al, th0, pv0
+  REAL,    DIMENSION(klon,klev) :: tpot_cen, tpot_edg, pvor_cen
+  REAL,    PARAMETER :: sg0=0.75  !--- Start level for PV=cte search loop
+  INTEGER, PARAMETER :: nadj=3    !--- Adjacent levs nb for thresholds detection
+  REAL,    PARAMETER :: w(5)=[0.1,0.25,0.3,0.25,0.1] !--- Vertical smoothing
+  INTEGER, SAVE :: k0
+  LOGICAL, SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(k0,first)
 !-------------------------------------------------------------------------------
   sub='dyn_tropopause'
@@ -42,5 +46,4 @@
   CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
   CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
-  CALL assert(SIZE(presnivs,1)==klev, TRIM(sub)//" presnivs klev")
   CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
   CALL assert(SHAPE(pplay)==[klon,klev  ],TRIM(sub)//" pplay shape")
@@ -48,57 +51,81 @@
 
   !--- DEFAULT THRESHOLDS
-  tpo0=380.;   IF(PRESENT(tpot0)) tpo0=tpot0  !--- In kelvins
-  pvo0=2.0;    IF(PRESENT(pvor0)) pvo0=pvor0  !--- In PVU
-  pmn0= 8000.; IF(PRESENT(pmin0)) pmn0=pmin0  !--- In pascals
-  pmx0=50000.; IF(PRESENT(pmax0)) pmx0=pmax0  !--- In pascals
-  kmin=klev-locate(presnivs(klev:1:-1),pmx0)+1
-  kmax=klev-locate(presnivs(klev:1:-1),pmn0)+1
+  th0=380.; IF(PRESENT(thet0)) th0=thet0   !--- In kelvins
+  pv0=  2.; IF(PRESENT(pvor0)) pv0=pvor0   !--- In PVU
+  IF(first) THEN
+    DO k0=1,klev; IF(aps(k0)/preff+bps(k0)<sg0) EXIT; END DO; first=.FALSE.
+  END IF
 
-  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS
-  DO k= 1, klev
-    DO i = 1, klon
-      t_pot(i,k) = t(i,k)*(preff/pplay(i,k))**RKAPPA
+  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES
+  DO i = 1,klon
+    tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA
+    tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA
+    DO k=2,klev
+      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
+      tpot_cen(i,k) =  t(i,k)                        *(preff/pplay(i,k))**RKAPPA
+      tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA
+      !--- FORCE QUANTITIES TO BE GROWING
+      IF(tpot_edg(i,k)<tpot_edg(i,k-1)) tpot_edg(i,k)=tpot_edg(i,k-1)+1.E-5
+      IF(tpot_cen(i,k)<tpot_cen(i,k-1)) tpot_cen(i,k)=tpot_cen(i,k-1)+1.E-5
     END DO
-  END DO
-
-  !--- TEMPERATURE AT CELLS INTERFACES (except in top layer)
-  t_edg(:,1) = ts(:)
-  DO k= 2, klev
-    DO i = 1, klon
-      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
-      t_edg(i,k) = t(i,k) + al*(t(i,k-1)-t(i,k))
-    END DO
+    !--- VERTICAL SMOOTHING
+    tpot_cen(i,:)=smooth(tpot_cen(i,:),w)
+    tpot_edg(i,:)=smooth(tpot_edg(i,:),w)
   END DO
 
   !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
-  DO k= 1, klev-1
-    DO i = 1, klon
-      IF(paprs(i,k)==paprs(i,k+1)) THEN; pvort(i,k)=HUGE(1.); CYCLE; END IF
-      pvort(i,k) = -1.E6 * RG * 2.*ROMEGA*ABS(SIN(latitude_deg(i)*RPI/180.))   &
-                         * ( t_edg(i,k+1)*(preff/paprs(i,k+1) )**RKAPPA        &
-                           - t_edg(i,k  )*(preff/paprs(i,k  ) )**RKAPPA)       &
-                         / ( paprs(i,k+1)  -     paprs(i,k  ) )
+  DO i = 1, klon
+    DO k= 1, klev-1
+      pvor_cen(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))&
+                   * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k))
     END DO
+    !--- VERTICAL SMOOTHING
+    pvor_cen(i,1:klev-1)=smooth(pvor_cen(i,1:klev-1),w)
   END DO
 
-  !--- LOCATE TROPOPAUSE
-  ptrop(:)=0.
+  !--- LOCATE TROPOPAUSE: LOWEST POINT BETWEEN THETA=380K AND PV=2PVU SURFACES.
   DO i = 1, klon
-    !--- Dynamical tropopause
-    it=kmax; DO WHILE(pvort(i,it)>pvo0.AND.it>=kmin); it=it-1; END DO
-    IF(pvort(i,it+1)/=pvort(i,it).AND.ALL([kmax,kmin-1]/=it) &
-      .AND.ALL([pvort(i,it),pvort(i,it+1)]/=HUGE(1.))) THEN
-      al = (pvort(i,it+1)-pvo0)/(pvort(i,it+1)-pvort(i,it))
-      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
-    END IF
-    !--- Potential temperature iso-surface
-    it = kmin-1+locate(t_pot(i,kmin:kmax),tpo0)
-    IF(t_pot(i,it+1)/=t_pot(i,it).AND.ALL([kmin-1,kmax]/=it)) THEN
-      al = (t_pot(i,it+1)-tpo0)/(t_pot(i,it+1)-t_pot(i,it))
-      ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
-    END IF
+    !--- UPPER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM TOP
+    DO kt=klev-1,1,-1; IF(ALL(ABS(pvor_cen(i,kt-nadj:kt))<=pv0)) EXIT; END DO
+    !--- LOWER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM BOTTOM
+    DO kb=k0,klev-1;   IF(ALL(ABS(pvor_cen(i,kb:kb+nadj))> pv0)) EXIT; END DO; kb=kb-1
+    !--- ISO-THETA POINT: THETA=380K       STARTING FROM TOP
+    DO kp=klev-1,1,-1; IF(ALL(ABS(tpot_cen(i,kp-nadj:kp))<=th0)) EXIT; END DO
+    !--- CHOOSE BETWEEN LOWER AND UPPER TROPOPAUSE
+    IF(2*COUNT(ABS(pvor_cen(i,kb:kt))>pv0)>kt-kb+1) kt=kb
+    !--- PV-DEFINED TROPOPAUSE
+    al = (ABS(pvor_cen(i,kt+1))-pv0)/ABS(pvor_cen(i,kt+1)-pvor_cen(i,kt))
+    dyn_tropopause(i) = pplay(i,kt+1)*(pplay(i,kt)/pplay(i,kt+1))**al
+    !--- THETA=380K IN THE TROPICAL REGION
+    al = (tpot_cen(i,kp+1)-th0)/(tpot_cen(i,kp+1)-tpot_cen(i,kp))
+    dyn_tropopause(i) = MAX( pplay(i,kp+1)*(pplay(i,kp)/pplay(i,kp+1))**al,    &
+                            dyn_tropopause(i) )
   END DO
 
-END SUBROUTINE dyn_tropopause
+END FUNCTION dyn_tropopause
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION smooth(v,w)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)         :: v(:), w(:)
+  REAL, DIMENSION(SIZE(v)) :: smooth
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: nv, nw, k, kb, ke, lb, le
+!-------------------------------------------------------------------------------
+  nv=SIZE(v); nw=(SIZE(w)-1)/2
+  DO k=1,nv
+    kb=MAX(k-nw,1 ); lb=MAX(2+nw   -k,1)
+    ke=MIN(k+nw,nv); le=MIN(1+nw+nv-k,1+2*nw)
+    smooth(k)=SUM(v(kb:ke)*w(lb:le))/SUM(w(lb:le))
+  END DO
+
+END FUNCTION smooth
+!
+!-------------------------------------------------------------------------------
 
 END MODULE tropopause_m
