Index: LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90	(revision 3069)
+++ LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90	(revision 3086)
@@ -16,7 +16,6 @@
 !    to the LMDZ vertical grid.
 !  * 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.
+!     - "lat_in":  input file latitudes.
+!     - "Ptrp_ou": target grid (LMDZ) tropopause pressure.
 !   so that the tropopause is shifted to the position of the one of LMDZ.
 !  Note that the ozone quantity conservation is not ensured.
@@ -35,6 +34,11 @@
 !    except that the latitudes are in ascending order in the input file.
 !  * We assume that all the input fields have the same coordinates.
-!  * The target vertical LMDZ grid is the grid of layer boundaries.
-!  * Regridding in pressure is conservative, second order.
+!  * The target vertical LMDZ grid is the grid of layer centers.
+!  * Regridding in pressure can be:
+!    - Ploc=='I': pressures provided at Interfaces    => conservative 2nd order
+!         OK for ozone coefficients regridding in Cariolle routines.
+!    - Ploc=='C': pressures provides at cells Centers => classical linear       
+!         OK for ozone vertical regridding, especially when torpopause
+!      adjustment is activated, to avoid "strairs shape effect" on profiles.
 !  * All the fields are regridded as a single multi-dimensional array, so it
 !    saves CPU time to call this procedure once for several NetCDF variables
@@ -43,6 +47,6 @@
 !    1) read from the file if "tropopause_air_pressure" field is available.
 !    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:
+!    3) computed using "tro3" and a fixed threshold otherwise, constant or
+!    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.
@@ -52,26 +56,28 @@
 !  * 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)".
+!  * Fields(i,j,k,l) are at longitude-latitude-name "rlonv(i)-rlatu(j)-nam(l)",
+!    pressure level/interval (Ploc=="C"/"I") "pcen_in(k)"/"pcen_in(k:k+1)]".
 !  * 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
+!  law in a limited zone made of 2 layers:
+!    1) between lower bound (lower than lowest tropopause) and LMDZ tropopause
+!    2) between LMDZ tropopause and upper bound (higher thzn highest tropopause)
 !  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.
+!        Sigma_str = Sigma^(1+alpha*phi(Sigma)), where:
+!   * alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
+!     This value shifts the file tropopause to the height of the one of LMDZ.
+!   * phi is quasi-linear (sections of 1/log function) in the adjacent intervals:
+!       - from 0 to 1 in [Sig_top,SigT_ou]
+!       - from 1 to 0 in [SigT_ou,Sig_bot]
+!  This quasi-triangular localization function ponderates alpha-law from one near
+!  the tropopause to zero each side apart.
 !
 ! * The following fields are on the global "dynamics" grid, as read from files:
-  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
-  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next ozone fields
-  REAL,    SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure
+  REAL,    SAVE, ALLOCATABLE :: v1 (:,:,:,:)       !--- Current  time ozone fields
+! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
+  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous time ozone fields
+  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next     time ozone fields
+  REAL,    SAVE, ALLOCATABLE :: pgm(:,:), pgp(:,:) !--- Ground     pressure
   REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
   REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
@@ -82,15 +88,15 @@
 !      * for monthly input files: julien is in [time_in(irec),time_in(irec+1)]
   LOGICAL, SAVE :: linterp                         !--- Interpolation in time
-  LOGICAL, SAVE :: lPrSurf                         !--- Surface pressure flag
-  LOGICAL, SAVE :: lPrTrop                         !--- Tropopause pressure flag
-  LOGICAL, SAVE :: lO3Trop                         !--- Tropopause ozone flag
+  LOGICAL, SAVE :: lPrSfile                        !--- Surface pressure flag
+  LOGICAL, SAVE :: lPrTfile                        !--- Tropopause pressure flag
+  LOGICAL, SAVE :: lO3Tfile                        !--- Tropopause ozone flag
   LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
 !$OMP THREADPRIVATE(lfirst)
   REAL,    PARAMETER :: pTropUp=9.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 parametrized O3 vmr at tropopause
-  LOGICAL, PARAMETER :: debug=.FALSE. !--- Force writefield_phy multiple outputs
+  REAL,    PARAMETER :: gamm   =0.4   !--- Max. stretched layer sigma ratio
+  REAL,    PARAMETER :: rho    =1.4   !--- Max tropopauses sigma ratio
+  REAL,    PARAMETER :: o3t0   =1.E-7 !--- Nominal O3 vmr at tropopause
+  LOGICAL, PARAMETER :: lO3Tpara=.FALSE. !--- Parametrized O3 vmr at tropopause
+  LOGICAL, PARAMETER :: ldebug=.TRUE.!--- Force writefield_phy multiple outputs
   REAL, PARAMETER :: ChemPTrMin=9.E+3 !--- Thresholds for minimum and maximum
   REAL, PARAMETER :: ChemPTrMax=4.E+4 !    chemical  tropopause pressure (Pa).
@@ -102,6 +108,6 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3,             &
-                             time_in, lon_in, lat_in, pcen_in, ptrop_ou)
+SUBROUTINE regr_pr_time_av(fID, nam, julien, Ploc, Pre_in, Pre_ou, v3, Pgnd_ou,&
+                                             time_in, lon_in, lat_in, Ptrp_ou)
 !
 !-------------------------------------------------------------------------------
@@ -118,22 +124,24 @@
   USE slopes_m,       ONLY: slopes
   USE mod_phys_lmdz_mpi_data,       ONLY: is_mpi_root
-  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev
+  USE mod_grid_phy_lmdz, ONLY: nlon=>nbp_lon, nlat=>nbp_lat, nlev_ou=>nbp_lev
   USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter
   USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
 !-------------------------------------------------------------------------------
 ! Arguments:
-  INTEGER, INTENT(IN)  :: fID                 !--- NetCDF file ID
+  INTEGER,           INTENT(IN) :: fID        !--- NetCDF file ID
   CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
-  REAL,    INTENT(IN)  :: julien              !--- Days since Jan 1st
-  !--- File & LMDZ (resp. decreasing & increasing order) interfaces pressure, Pa
-  REAL,    INTENT(IN)  :: pint_in(:)          !--- in:  file          (nlev_in+1)
-  REAL,    INTENT(IN)  :: pint_ou(:,:)        !--- out: LMDZ     (klon,nbp_lev+1)
-  REAL,    INTENT(OUT) :: v3(:,:,:)           !--- Regridded fld (klon,llm,n_var)
+  REAL,              INTENT(IN) :: julien     !--- Days since Jan 1st
+  CHARACTER(LEN=1),  INTENT(IN) :: Ploc       !--- Pressures locations
+  !--- File/LMDZ (resp. decreasing & increasing order) pressure, Pa
+  !    At cells centers or interfaces depending on "Ploc" keyword (C/I)
+  REAL,    INTENT(IN)  :: Pre_in(:)           !--- in:  file      (nlev_in[+1])
+  REAL,    INTENT(IN)  :: Pre_ou(:,:)         !--- out: LMDZ (klon,nlev_ou[+1])
+  REAL,    INTENT(OUT) :: v3(:,:,:)           !--- Regr. fld (klon,nlev_ou,n_var)
+  REAL,    INTENT(IN), OPTIONAL :: Pgnd_ou(:) !--- LMDZ ground pressure   (klon)
   REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
                                               !    since Jan 1 of current year
   REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- File longitudes vector (klon)
   REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- File latitudes  vector (klon)
-  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- File mid-layers pres (nlev_in)
-  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- LMDZ tropopause pres   (klon)
+  REAL,    INTENT(IN), OPTIONAL :: Ptrp_ou(:) !--- LMDZ tropopause pres   (klon)
 !-------------------------------------------------------------------------------
 ! Local variables:
@@ -142,45 +150,56 @@
   CHARACTER(LEN=80)  :: sub
   CHARACTER(LEN=320) :: str
-  INTEGER :: vID, ncerr, n_var, ibot, ibo0, nn, itrp
-  INTEGER :: i, nlev_in, n_dim, itop, ito0, i0
+  INTEGER :: vID, ncerr, n_var, ibot, iout, nn
+  INTEGER :: i, nlev_in, n_dim, itop, itrp, i0
   LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
   REAL    :: y_frac                           !--- Elapsed year fraction
   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))
-! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
-! v2: Field scattered to the partial "physics" horizontal grid.
+  REAL    :: Sig_bot, Sig_top                 !--- Bounds of quasi-hat function
+  REAL    :: Sig_bo0, Sig_to0                 !--- Lower/upper tropopauses
+  REAL    :: Sig_in (SIZE(Pre_in))            !--- Input field sigma levels
+  REAL    :: Sig_ou (SIZE(Pre_ou,2))          !--- Output LMDZ sigma levels
+  REAL    :: phi    (SIZE(Pre_ou,2))          !--- Stretching exponent anomaly
+  REAL    :: Pstr_ou(SIZE(Pre_ou,2))          !--- Stretched pressure levels
+  REAL    :: Pres_ou(SIZE(Pre_ou,2))          !--- Pre_ou(i,:), reversed order
+  REAL, DIMENSION(nlon, nlat) :: pg1,      pt1,      ot1
+  REAL, DIMENSION(klon)       :: Pgnd_in,  Ptrp_in,  Otrp_in
+  REAL, DIMENSION(klon)       :: Ptrop_ou, Pgrnd_ou
+! * The following fields are scattered to the partial "physics" horizontal grid.
+  REAL, POINTER :: v2(:,:,:)                  !--- Current  time ozone fields
 !     In the 2D file case, the values are the same for all longitudes.
-!     "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, ptropou
+!     "v2(i, k, l)" is at longitude-latitude "xlon(i)-xlat(i)" and name "nam(l)"
+! Both are:          * if Ploc=='I' in pressure interval "press_in_edg(k:k+1)"
+!                    * if Ploc=='C' at pressure          "press_in_cen(k)"
+  REAL, TARGET :: &
+    v2i(klon,SIZE(Pre_in)-1,SIZE(nam)), &     !--- v2 in Ploc=='I' case
+    v2c(klon,SIZE(Pre_in)  ,SIZE(nam))        !--- v2 in Ploc=='C' case
   LOGICAL :: ll
 !--- For debug
-  REAL, DIMENSION(klon)           :: ptrop_in, ptrop_ef!, dzStrain, dzStrain0
-! REAL, DIMENSION(klon,nbp_lev+1) :: pstrn_ou, phii
+  REAL, DIMENSION(klon)             :: Ptrop_in, Ptrop_ef
+  REAL, DIMENSION(klon)             :: dzStrain, dzStrain0
+  REAL, DIMENSION(klon,SIZE(Pre_ou,2)) :: Pstrn_ou, phii
 !-------------------------------------------------------------------------------
   sub="regr_pr_time_av"
-  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")
+  nlev_in=SIZE(Pre_in); IF(Ploc=='I') nlev_in=nlev_in-1
+  IF(Ploc=='I') THEN; v2 => v2i; ELSE; v2 => v2c; END IF
+  CALL assert(SIZE(Pre_ou,1)==klon,TRIM(sub)//" Pre_ou klon")
+  CALL assert(SIZE(v3,1)==klon,    TRIM(sub)//" v3 klon")
+  CALL assert(SIZE(v3,2)==nlev_ou, TRIM(sub)//" v3 nlev_ou")
+  IF(Ploc=='I') CALL assert(SIZE(Pre_ou,2)==nlev_ou+1,TRIM(sub)//" Pre_ou nlev_ou+1")
+  IF(Ploc=='C') CALL assert(SIZE(Pre_ou,2)==nlev_ou  ,TRIM(sub)//" Pre_ou nlev_ou")
   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")
-  lAdjTro=PRESENT(ptrop_ou)
-  IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) &
-  CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1)
+  IF(PRESENT(Pgnd_ou)) CALL assert(SIZE(Pgnd_ou)==klon,TRIM(sub)//" Pgnd_ou klon")
+  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(Ptrp_ou)) CALL assert(SIZE(Ptrp_ou)==klon,TRIM(sub)//" Ptrp_ou klon")
+  lAdjTro=PRESENT(Ptrp_ou)
+  IF(lAdjTro) THEN
+    IF(.NOT.PRESENT(lat_in)) &
+      CALL abort_physic(sub, 'Missing lat_in (required if adjust_tropopause=T)', 1)
+    IF(.NOT.PRESENT(Pgnd_ou).AND.Ploc=='C') &
+      CALL abort_physic(sub, 'Missing ground Pr(required if adjust_tropopause=T)', 1)
+    IF(PRESENT(Pgnd_ou)) THEN; Pgrnd_ou=Pgnd_ou; ELSE; Pgrnd_ou=Pre_ou(:,1); END IF
+  END IF
 
   !$OMP MASTER
@@ -189,31 +208,27 @@
     !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE
     IF(lfirst) THEN
-      lPrSurf=NF90_INQ_VARID(fID,"ps"                     ,vID)==NF90_NOERR
-      lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
-      lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
+      lPrSfile=lAdjTro.AND.NF90_INQ_VARID(fID,"ps"                     ,vID)==NF90_NOERR
+      lPrTfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
+      lO3Tfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
       CALL NF95_INQ_DIMID(fID,"time",vID)
       CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
       linterp=PRESENT(time_in).AND.ntim_in==14
+      ALLOCATE(v1(nlon,nlat,nlev_in,n_var))
       IF(linterp) THEN
-        ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
-        ALLOCATE(v1p(nbp_lon,nbp_lat,nlev_in,n_var))
-        ALLOCATE(psm(nbp_lon,nbp_lat),psp(nbp_lon,nbp_lat))
-        ALLOCATE(ptm(nbp_lon,nbp_lat),ptp(nbp_lon,nbp_lat))
-        IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat))
+        ALLOCATE(v1m(nlon,nlat,nlev_in,n_var),v1p(nlon,nlat,nlev_in,n_var))
+        IF(lPrSfile) ALLOCATE(pgm(nlon,nlat),pgp(nlon,nlat))
+        IF(lPrTfile) ALLOCATE(ptm(nlon,nlat),ptp(nlon,nlat))
+        IF(lO3Tfile) ALLOCATE(otm(nlon,nlat),otp(nlon,nlat))
       END IF
       !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
-      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 IF(lo3tp) THEN
-        CALL msg(lAdjTro,'    PARAMETRIZED O3 concentration at tropopause.')
-      ELSE
-        CALL msg(lAdjTro,'    CONSTANT O3 concentration at tropopause.')
-      END IF
+      IF(lAdjTro) itrp0=locate(Pre_in,pTropUp)
+      CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.'    ,sub)
+      CALL msg(lPrSfile,'Using GROUND PRESSURE from input O3 forcing file.',sub)
+      CALL msg(lAdjTro ,'o3 forcing file tropopause location uses:'        ,sub)
+      IF(lPrTfile)      THEN; str='    INPUT FILE PRESSURE'
+      ELSE IF(lO3Tfile) THEN; str='    INPUT FILE O3 CONCENTRATION'
+      ELSE IF(lO3Tpara) THEN; str='    PARAMETRIZED O3 concentration'
+      ELSE;                   str='    CONSTANT O3 concentration'; END IF
+      CALL msg(lAdjTro,TRIM(str)//' at tropopause')
     END IF
 
@@ -223,35 +238,41 @@
     !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
     IF(linterp) THEN
-      WRITE(str,'(3(a,f12.8))')'Interpolating O3 at julian day ',julien,' from '&
-      &//'fields at times ',time_in(irec),' and ', time_in(irec+1)
+      WRITE(str,'(a,f12.8,2(a,f5.1))')'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
-      IF(lPrSurf) ps1=al*psm+(1.-al)*psp
-      IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp
-      IF(lO3Trop) ot1=al*otm+(1.-al)*otp
+      IF(lPrSfile) pg1=al*pgm+(1.-al)*pgp
+      IF(lPrTfile) pt1=al*ptm+(1.-al)*ptp
+      IF(lO3Tfile) ot1=al*otm+(1.-al)*otp
     END IF
   END IF
   !$OMP END MASTER
   IF(lfirst) THEN
-    lfirst=.FALSE.;      CALL bcast(lfirst)
-    IF(lAdjTro)          CALL bcast(itrp0)
-    CALL bcast(lPrTrop); CALL bcast(lPrSurf)
-    CALL bcast(lO3Trop); CALL bcast(linterp)
+    lfirst=.FALSE.;       CALL bcast(lfirst)
+    IF(lAdjTro)           CALL bcast(itrp0)
+    CALL bcast(lPrSfile); CALL bcast(lPrTfile)
+    CALL bcast(lO3Tfile); 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(lPrTrop) CALL scatter2d(pt1,pt2)
-  IF(lO3Trop) CALL scatter2d(ot1,ot2)
-
-  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
-  IF(.NOT.lAdjTro) THEN
+  IF(lPrSfile) CALL scatter2d(pg1,Pgnd_in)
+  IF(lPrTfile) CALL scatter2d(pt1,Ptrp_in)
+  IF(lO3Tfile) CALL scatter2d(ot1,Otrp_in)
+  !--- No ground pressure in input file => choose it to be the one of LMDZ
+  IF(lAdjTro.AND..NOT.lPrSfile) Pgnd_in(:)=Pgrnd_ou(:)
+  
+!-------------------------------------------------------------------------------
+  IF(.NOT.lAdjTro) THEN       !--- REGRID IN PRESSURE ; NO TROPOPAUSE ADJUSTMENT
+!-------------------------------------------------------------------------------
     DO i=1,klon
-      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(:)))
+      Pres_ou=Pre_ou(i,SIZE(Pre_ou,2):1:-1)   !--- pplay & paprs are decreasing
+      IF(Ploc=='C') CALL regr_lint   (1,v2(i,:,:), LOG(Pre_in(:)),             &
+        LOG(Pres_ou(:)), v3(i,nlev_ou:1:-1,:))
+      IF(Ploc=='I') CALL regr_conserv(1,v2(i,:,:),     Pre_in(:) ,             &
+            Pres_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
     END DO
-  ELSE
+!-------------------------------------------------------------------------------
+  ELSE                        !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT
+!-------------------------------------------------------------------------------
     y_frac=(REAL(days_elapsed)+jH_cur)/year_len
 
@@ -259,8 +280,8 @@
     DO i=1,klon
 
-      !--- 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/OUTPUT (FILE/LMDZ) SIGMA LEVELS IN CURRENT COLUMN
+      Pres_ou   = Pre_ou(i,SIZE(Pre_ou,2):1:-1)!--- pplay & paprs are decreasing
+      Sig_in(:) = Pre_in (:)/Pgnd_in(i)            !--- increasing values
+      Sig_ou(:) = Pres_ou(:)/Pgnd_ou(i)            !--- increasing values
 
       !--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered
@@ -268,89 +289,90 @@
       ! ozone hole fooling the crude chemical tropopause detection algorithm.
       SigT_in = get_SigTrop(i,itrp)
-      SigT_in=MIN(SigT_in,ChemPTrMax/ps2(i))      !--- too low  value filtered
-      SigT_in=MAX(SigT_in,ChemPTrMin/ps2(i))      !--- too high value filtered
+      SigT_in=MIN(SigT_in,ChemPTrMax/Pgnd_in(i))   !--- too low  value filtered
+      SigT_in=MAX(SigT_in,ChemPTrMin/Pgnd_ou(i))   !--- too high value filtered
 
       !--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
-      ! dynamical tropopause (especially in filaments) are conterbalanced with a
-      ! filter ensuring it stays within a certain distance around input (file)
+      ! dynamical tropopause (especially in filaments) are conterbalanced with
+      ! a filter ensuring it stays within a certain distance around input (file)
       ! tropopause, hence avoiding avoid a too thick stretched region ; a final
       ! extra-safety filter keeps the tropopause pressure value realistic.
-      SigT_ou = ptrop_ou(i)/ps2(i)
-      IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho !--- too low  value w/r input
-      IF(SigT_ou>SigT_in*rho) SigT_ou=SigT_in*rho !--- too high value w/r input
-      SigT_ou=MIN(SigT_ou,DynPTrMax/ps2(i))       !--- too low  value filtered
-      SigT_ou=MAX(SigT_ou,DynPTrMin/ps2(i))       !--- too high value filtered
-      ptropou(i)=SigT_ou*ps2(i)
-
-      !--- FULLY STRETCHED LAYER BOUNDS (alpha power law fully applied)
+      SigT_ou = Ptrp_ou(i)/Pgnd_ou(i)
+      IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho  !--- too low  value w/r input
+      IF(SigT_ou>SigT_in*rho) SigT_ou=SigT_in*rho  !--- too high value w/r input
+      SigT_ou=MIN(SigT_ou,DynPTrMax/Pgnd_ou(i))    !--- too low  value filtered
+      SigT_ou=MAX(SigT_ou,DynPTrMin/Pgnd_ou(i))    !--- too high value filtered
+      Ptrop_ou(i)=SigT_ou*Pgnd_ou(i)
+      iout = locate(Sig_ou(:),SigT_ou)
+
+      !--- POWER LAW COEFFICIENT FOR TROPOPAUSES MATCHING
       alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
-      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 (fraction of alpha applied)
-      ! The used formulae ensure the first derivative stays positive.
-      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)
-
-      !--- STRETCHING POWER LAW FUNCTION:
-      !    0    in [0,Sig_to0] and [Sig_bo0,1] ;    1    in [Sig_bot,Sig_top]
-      !    0->1 in [Sig_to0,Sig_top]                1->0 in [Sig_bot,Sig_bo0]
+
+      !--- DETERMINE STRETCHING DOMAIN UPPER AND LOWER BOUNDS
+      Sig_bo0 = MAX(SigT_in,SigT_ou)               !--- lowest  tropopause
+      Sig_to0 = MIN(SigT_in,SigT_ou)               !--- highest tropopause
+      beta    = (Sig_bo0/Sig_to0)**gamm            !--- stretching exponent
+      Sig_bot = MIN(Sig_bo0*beta,0.1*(9.+Sig_bot)) !--- must be <1
+      ibot = locate(Sig_ou(:),Sig_bot)             !--- layer index
+      IF(ibot-iout<2) THEN                         !--- at least one layer thick
+        ibot=MIN(iout+2,nlev_ou); Sig_bot=Sig_ou(ibot)
+      END IF
+      Sig_top = Sig_to0/beta                       !--- upper bound
+      itop = locate(Sig_ou(:),Sig_top)             !--- layer index
+      IF(iout-itop<2) THEN                         !--- at least one layer thick
+        itop=MAX(iout-2,1); Sig_top=Sig_ou(itop)
+      END IF
+
+      !--- STRETCHING POWER LAW LOCALIZATION FUNCTION:
+      !    0 in [0,Sig_top]    0->1 in [Sig_top,SigT_ou]
+      !    0 in [Sig_bot,1]    1->0 in [SigT_ou, Sig_bot]
       phi(:)=0.
-      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)
+      phi(itop+1:iout) = (1.-LOG(Sig_top)/LOG(Sig_ou(itop+1:iout)))&
+                            *LOG(SigT_ou)/LOG(SigT_ou/Sig_top)
+      phi(iout+1:ibot) = (1.-LOG(Sig_bot)/LOG(Sig_ou(iout+1:ibot)))&
+                            *LOG(SigT_ou)/LOG(SigT_ou/Sig_bot)
 
       !--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
-      pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
+      Pstr_ou(:) = Pres_ou(:) * 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(:)))
+      IF(Ploc=='C') CALL regr_lint   (1, v2(i,:,:), LOG(Pre_in(:)),            &
+        LOG(Pstr_ou(:)), v3(i,nlev_ou:1:-1,:))
+      IF(Ploc=='I') CALL regr_conserv(1, v2(i,:,:),     Pre_in(:) ,            &
+            Pstr_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_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',  &
+      i0=nlev_ou-locate(Pres_ou(:),Ptrop_ou(i))+1
+      ll=check_ozone(v3(i, 1:i0-1   ,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', &
+      ll=check_ozone(v3(i,i0:nlev_ou,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)
 
-      IF(debug) THEN
-!        dzStrain0(i)=7.*LOG(Sig_bo0/Sig_to0)
-!        dzStrain (i)=7.*LOG(Sig_bot/Sig_top)
-        ptrop_in(i) = SigT_in*ps2(i)
-!        Pstrn_ou(i,:)= pstr_ou
-!        phii(i,:)=phi(:)
-        ptrop_ef(i)=chem_tropopause(i, itrp, locate(pintou(:),PTropUp),        &
-                           pintou(:), v3(:,nbp_lev+1:1:-1,1),o3trop=o3t0)
+      IF(ldebug) THEN
+        dzStrain0(i) = SIGN(7.*LOG(Sig_bo0/Sig_to0),SigT_in-SigT_ou)
+        dzStrain (i) = SIGN(7.*LOG(Sig_bot/Sig_top),SigT_in-SigT_ou)
+        Ptrop_in (i) = SigT_in*Pgnd_in(i)
+        Pstrn_ou(i,:)= Pstr_ou
+        phii(i,:)    = phi(:)
+        Ptrop_ef(i)  = PTrop_chem(i, itrp, locate(Pres_ou(:),PTropUp),    &
+                             Pres_ou(:), v3(:,nlev_ou:1:-1,1),o3trop=o3t0)
       END IF
     END DO
-    IF(debug) THEN
-!      CALL writefield_phy('PreSt_ou' ,Pstrn_ou,nbp_lev+1) !--- Strained pressure
-!      CALL writefield_phy('dzStrain' ,dzStrain ,1)   !--- Fully & total strained
-!      CALL writefield_phy('dzStrain0',dzStrain0,1)   !--- regions thickness
-!      CALL writefield_phy('phi',phii,nbp_lev+1)      !--- Localisation function
-      !--- Tropopauses pressures:
-      CALL writefield_phy('PreTr_in',ptrop_in,1)     !--- Input and effective
-      CALL writefield_phy('PreTr_ef',ptrop_ef,1)     !    chemical tropopauses
-    END IF
-  END IF
-  IF(debug) THEN
-    IF(lAdjTro) &
-    CALL writefield_phy('PreTr_ou',ptropou,1)        !--- LMDz dyn tropopause
+  END IF
+  IF(ldebug.AND.lAdjTro) THEN
+    CALL writefield_phy('PreSt_ou' ,Pstrn_ou,SIZE(Pre_ou,2)) !--- Strained Pres
+    CALL writefield_phy('dzStrain' ,dzStrain ,1)     !--- Strained thickness
+    CALL writefield_phy('dzStrain0',dzStrain0,1)     !--- Tropopauses distance
+    CALL writefield_phy('phi',phii,nlev_ou)          !--- Localization function
+    !--- Tropopauses pressures:
+    CALL writefield_phy('PreTr_in',Ptrop_in,1)       !--- Input and effective
+    CALL writefield_phy('PreTr_ou',Ptrop_ou,1)       !--- LMDz dyn tropopause
+    CALL writefield_phy('PreTr_ef',Ptrop_ef,1)       !--- Effective chem tropop
+  END IF
+  IF(ldebug) THEN
     CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
-    CALL writefield_phy('Ozone_ou',v3(:,:,1),nbp_lev)!--- Output ozone field
-    CALL writefield_phy('Pint_ou' ,pint_ou,1+nbp_lev)!--- LMDZ unstreched Press
+    CALL writefield_phy('Ozone_ou',v3(:,:,1),nlev_ou)!--- Output ozone field
+    CALL writefield_phy('Pres_ou' ,Pre_ou,SIZE(Pre_ou,2))!--- LMDZ Pressure
   END IF
 
@@ -369,7 +391,7 @@
     CALL get_3Dfields(v1)               !--- Read ozone field(s)
     IF(lAdjTro) THEN                    !--- Additional files for fields strain
-      IF(lPrSurf) CALL get_2Dfield(ps1,"ps")
-      IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure")
-      IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause")
+      IF(lPrSfile) CALL get_2Dfield(pg1,"ps")
+      IF(lPrTfile) CALL get_2Dfield(pt1,"tropopause_air_pressure")
+      IF(lO3Tfile) CALL get_2Dfield(ot1,"tro3_at_tropopause")
     END IF
   ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
@@ -383,7 +405,7 @@
       CALL get_3Dfields(v1m)            !--- Read ozone field(s)
       IF(lAdjTro) THEN                  !--- Additional files for fields strain
-        IF(lPrSurf) CALL get_2Dfield(psm,"ps")
-        IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure")
-        IF(lO3Trop) CALL get_2Dfield(otm,"tro3_at_tropopause")
+        IF(lPrSfile) CALL get_2Dfield(pgm,"ps")
+        IF(lPrTfile) CALL get_2Dfield(ptm,"tropopause_air_pressure")
+        IF(lO3Tfile) CALL get_2Dfield(otm,"tro3_at_tropopause")
       END IF
     ELSE                                !=== SHIFT FIELDS
@@ -394,7 +416,7 @@
       v1m=v1p                           !--- Ozone fields
       IF(lAdjTro) THEN                  !--- Additional files for fields strain
-        IF(lPrSurf) psm=psp             !--- Surface pressure
-        IF(lPrTrop) ptm=ptp             !--- Tropopause pressure
-        IF(lO3Trop) otm=otp             !--- Tropopause ozone
+        IF(lPrSfile) pgm=pgp             !--- Surface pressure
+        IF(lPrTfile) ptm=ptp             !--- Tropopause pressure
+        IF(lO3Tfile) otm=otp             !--- Tropopause ozone
       END IF
     END IF
@@ -405,7 +427,7 @@
     CALL get_3Dfields(v1p)              !--- Read ozone field(s)
     IF(lAdjTro) THEN                    !--- Additional files for fields strain
-      IF(lPrSurf) CALL get_2Dfield(psp,"ps")
-      IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure")
-      IF(lO3Trop) CALL get_2Dfield(otp,"tro3_at_tropopause")
+      IF(lPrSfile) CALL get_2Dfield(pgp,"ps")
+      IF(lPrTfile) CALL get_2Dfield(ptp,"tropopause_air_pressure")
+      IF(lO3Tfile) CALL get_2Dfield(otp,"tro3_at_tropopause")
     END IF
     irec=irec-1
@@ -438,8 +460,8 @@
   !--- Flip latitudes: ascending in input file, descending in "rlatu".
   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
+    v(1,:) = v(1,nlat:1:-1)
+    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nlon-1)  !--- Duplication
   ELSE
-    v(:,:) = v(:,nbp_lat:1:-1)
+    v(:,:) = v(:,nlat:1:-1)
   END IF
 
@@ -471,8 +493,8 @@
   !--- Flip latitudes: ascending in input file, descending in "rlatu".
   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
+    v(1,:,:,:) = v(1,nlat:1:-1,:,:)
+    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1)  !--- Duplication
   ELSE
-    v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:)
+    v(:,:,:,:) = v(:,nlat:1:-1,:,:)
   END IF
 
@@ -485,25 +507,20 @@
 !-------------------------------------------------------------------------------
 !
-FUNCTION get_SigTrop(ih,it)
-!
-!-------------------------------------------------------------------------------
-! Arguments:
+FUNCTION get_SigTrop(ih,it) RESULT(out)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL                 :: out
   INTEGER, INTENT(IN)  :: ih
   INTEGER, INTENT(OUT) :: it
-  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)
+!-------------------------------------------------------------------------------
+  !--- Pressure at tropopause read from the forcing file
+       IF(lPrTfile) THEN; out=Ptrp_in(ih)/Pgnd_in(ih); RETURN; END IF
+
+  !--- Chemical tropopause definition based on a particular threshold
+       IF(lO3Tfile) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),Otrp_in(ih))
+  ELSE IF(lO3Tpara) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1))
+  ELSE                  ; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),o3t0); END IF
+  out=out/Pgnd_in(ih)
 
 END FUNCTION get_SigTrop
@@ -514,5 +531,5 @@
 !-------------------------------------------------------------------------------
 !
-FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop)
+FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out)
 !
 !-------------------------------------------------------------------------------
@@ -529,17 +546,16 @@
 !-------------------------------------------------------------------------------
 ! Arguments:
-  REAL ::                   chem_tropopause    !--- Pressure at tropopause
+  REAL                        :: out           !--- 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)  :: pres(:)       !--- Pressure profile, 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:
-  REAL    :: o3t                               !--- Ozone concent. at tropopause
-  REAL    :: al                                !--- Interpolation coefficient
-  REAL    :: coef                              !--- Coeff of latitude modulation
+  REAL :: o3t                                  !--- Ozone concent. at tropopause
+  REAL :: al                                   !--- Interpolation coefficient
+  REAL :: coef                                 !--- Coeff of latitude modulation
   REAL, PARAMETER :: co3(3)=[91.,28.,20.]      !--- Coeff for o3 at tropopause
 !-------------------------------------------------------------------------------
@@ -556,19 +572,16 @@
   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)
+  IF(Ploc=='C') out =      pres(it)**(1.-al) * pres(it+1)**al
+  IF(Ploc=='I') out = SQRT(pres(it)**(1.-al) * pres(it+2)**al *pres(it+1))
+  it = locate(pres(:), out)                    !--- pres(it)<Ptrop<pres(it+1)
+
+END FUNCTION PTrop_chem
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out)
 !
 !-------------------------------------------------------------------------------
@@ -576,5 +589,5 @@
 !-------------------------------------------------------------------------------
 ! Arguments:
-  LOGICAL                      :: check_ozone      !--- .T. => some wrong values
+  LOGICAL                      :: out          !--- .T. => some wrong values
   REAL,             INTENT(IN) :: o3col(:), lon, lat
   INTEGER,          INTENT(IN) :: ilev0
@@ -592,6 +605,5 @@
   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
-  IF(prt_level<1) RETURN
+  out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN
 
   !--- SOME TOO LOW VALUES FOUND
