Changeset 3086


Ignore:
Timestamp:
Nov 22, 2017, 11:01:27 PM (6 years ago)
Author:
dcugnet
Message:

Modification of the ozone interpolation routine to ensure more realistic
profiles, in particular when the tropopause adjustment is activated, as
well as an improved accuracy on the chemical tropopause of streched ozone
field:

  • stretching function has no longer a trapezoidal shape, but is a hat

function, which top is at the target (ldmz) tropopause, with lower /
upper bounds respectively lower / higher than the lowest / highest
tropopause (lmdz and file). This helps to avoid stretching the profile
too far from the strictly required zone.

  • revert to a classical linear vertical interpolation on log-pressure

because the stairs-shape effect of the conservative interpolation can not
be avoided, even at second order method. Note that the second order
conservative method is still used for Cariolle linearized ozone
coefficients interpolation: the regr_pr_time routine now switches from
one method to the other depending on flag "Ploc" (pressures must be given
at layers interfaces - Ploc=='I' - for conservative method and layers
centers - Ploc=='C' - for the regular one).

  • force the stretching function to be discretized on at least three

levels to avoid the effect of the stretching to be reduced to a single
point shifting of the ozone profile when discretized.

Location:
LMDZ6/trunk/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r3082 r3086  
    20712071          IF(adjust_tropopause) THEN
    20722072             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
    2073                       ro3i,  press_edg_climoz,  paprs,   wo, time_climoz,    &
    2074                       longitude_deg,   latitude_deg,    press_cen_climoz,    &
     2073                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
     2074                      time_climoz ,  longitude_deg,   latitude_deg,          &
    20752075                      dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
    20762076          ELSE
    2077              CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
    2078                       ro3i,  press_edg_climoz,  paprs,  wo,  time_climoz)
     2077             CALL regr_pr_time_av(ncid_climoz,  vars_climoz(1:read_climoz),  &
     2078                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
     2079                      time_climoz )
    20792080          END IF
    20802081          ! Convert from mole fraction of ozone to column density of ozone in a

  • LMDZ6/trunk/libf/phylmd/regr_pr_comb_coefoz_m.F90

    r3043 r3086  
    129129
    130130    call regr_pr_time_av(ncid, (/"a2           ", "a4           ", "a6           ", &
    131          "P_net_Mob    ", "r_Mob        ", "temp_Mob     ", "R_Het        "/), REAL(julien-1), &
    132          press_in_edg, paprs, coefoz)
     131         "P_net_Mob    ", "r_Mob        ", "temp_Mob     ", "R_Het        "/),     &
     132         REAL(julien-1), 'I', press_in_edg, paprs, coefoz)
    133133    a2 = coefoz(:, :, 1)
    134134    a4_mass = coefoz(:, :, 2) * 48. / 29.

  • LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90

    r3069 r3086  
    1616!    to the LMDZ vertical grid.
    1717!  * the forcing fields are stretched if the following arguments are present:
    18 !     - "lat_in":   input file latitudes.
    19 !     - "pcen_in":  input file cells center pressure levels.
    20 !     - "ptrop_ou": target grid (LMDZ) tropopause pressure.
     18!     - "lat_in":  input file latitudes.
     19!     - "Ptrp_ou": target grid (LMDZ) tropopause pressure.
    2120!   so that the tropopause is shifted to the position of the one of LMDZ.
    2221!  Note that the ozone quantity conservation is not ensured.
    … …  
    3534!    except that the latitudes are in ascending order in the input file.
    3635!  * We assume that all the input fields have the same coordinates.
    37 !  * The target vertical LMDZ grid is the grid of layer boundaries.
    38 !  * Regridding in pressure is conservative, second order.
     36!  * The target vertical LMDZ grid is the grid of layer centers.
     37!  * Regridding in pressure can be:
     38!    - Ploc=='I': pressures provided at Interfaces    => conservative 2nd order
     39!         OK for ozone coefficients regridding in Cariolle routines.
     40!    - Ploc=='C': pressures provides at cells Centers => classical linear       
     41!         OK for ozone vertical regridding, especially when torpopause
     42!      adjustment is activated, to avoid "strairs shape effect" on profiles.
    3943!  * All the fields are regridded as a single multi-dimensional array, so it
    4044!    saves CPU time to call this procedure once for several NetCDF variables
    … …  
    4347!    1) read from the file if "tropopause_air_pressure" field is available.
    4448!    2) computed using "tro3" and "tro3_at_tropopause' (if available).
    45 !    3) computed using "tro3" and a fixed threshold otherwise, determined using
    46 !    an empirical three parameters law:
     49!    3) computed using "tro3" and a fixed threshold otherwise, constant or
     50!    determined using an empirical three parameters law:
    4751!         o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3)
    4852!       => co1 and co2 are in ppbV, and co3 in degrees.
    … …  
    5256!  * Fields with suffix "m"/"p" are at the closest records earlier/later than
    5357!  the mid-julian day "julien", on the global "dynamics" horizontal grid.
    54 !  * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
    55 !    interval "pint_in(k:k+1)]" and variable "nam(l)".
     58!  * Fields(i,j,k,l) are at longitude-latitude-name "rlonv(i)-rlatu(j)-nam(l)",
     59!    pressure level/interval (Ploc=="C"/"I") "pcen_in(k)"/"pcen_in(k:k+1)]".
    5660!  * In the 2D file case, the values are the same for all longitudes.
    5761!  * The tropopause correction works like this: the input fields (file) are
    5862!  interpolated on output (LMDZ) pressure field, which is streched using a power
    59 !  law in a limited zone made of 3 layers:
    60 !    1) between the two tropopauses (file and LMDZ)
    61 !    2) in an upper and a lower transitions layers
     63!  law in a limited zone made of 2 layers:
     64!    1) between lower bound (lower than lowest tropopause) and LMDZ tropopause
     65!    2) between LMDZ tropopause and upper bound (higher thzn highest tropopause)
    6266!  The stretching function has the following form:
    63 !    Sigma_str = Sigma^(1+alpha), with alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
    64 !  This value shifts the file tropopause to the height of the one of LMDZ.
    65 !  The stretching is fully applied in the central zone only, and only partially
    66 !  in the transitions zones, thick enough to guarantee a growing stretched
    67 !  pressure field. The ponderation function for alpha to modulate the stretching
    68 !  is constant equal to 1 in the central layer, and quasi-linear (from 1 to 0)
    69 !  in the transition layers (from 1 to 0 ; in fact: sections of 1/log function),
    70 !  making this localisation function quasi-trapezoidal.
     67!        Sigma_str = Sigma^(1+alpha*phi(Sigma)), where:
     68!   * alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
     69!     This value shifts the file tropopause to the height of the one of LMDZ.
     70!   * phi is quasi-linear (sections of 1/log function) in the adjacent intervals:
     71!       - from 0 to 1 in [Sig_top,SigT_ou]
     72!       - from 1 to 0 in [SigT_ou,Sig_bot]
     73!  This quasi-triangular localization function ponderates alpha-law from one near
     74!  the tropopause to zero each side apart.
    7175!
    7276! * The following fields are on the global "dynamics" grid, as read from files:
    73   REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
    74   REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next ozone fields
    75   REAL,    SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure
     77  REAL,    SAVE, ALLOCATABLE :: v1 (:,:,:,:)       !--- Current  time ozone fields
     78! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
     79  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous time ozone fields
     80  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next     time ozone fields
     81  REAL,    SAVE, ALLOCATABLE :: pgm(:,:), pgp(:,:) !--- Ground     pressure
    7682  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
    7783  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
    … …  
    8288!      * for monthly input files: julien is in [time_in(irec),time_in(irec+1)]
    8389  LOGICAL, SAVE :: linterp                         !--- Interpolation in time
    84   LOGICAL, SAVE :: lPrSurf                         !--- Surface pressure flag
    85   LOGICAL, SAVE :: lPrTrop                         !--- Tropopause pressure flag
    86   LOGICAL, SAVE :: lO3Trop                         !--- Tropopause ozone flag
     90  LOGICAL, SAVE :: lPrSfile                        !--- Surface pressure flag
     91  LOGICAL, SAVE :: lPrTfile                        !--- Tropopause pressure flag
     92  LOGICAL, SAVE :: lO3Tfile                        !--- Tropopause ozone flag
    8793  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
    8894!$OMP THREADPRIVATE(lfirst)
    8995  REAL,    PARAMETER :: pTropUp=9.E+3 !--- Value  <  tropopause pressure (Pa)
    90   REAL,    PARAMETER :: gamm = 0.4    !--- Relative thickness of transitions
    91   REAL,    PARAMETER :: rho  = 1.3    !--- Max tropopauses sigma ratio
    92   REAL,    PARAMETER :: o3t0 = 1.E-7 !--- Nominal O3 vmr at tropopause
    93   LOGICAL, PARAMETER :: lo3tp=.FALSE. !--- Use parametrized O3 vmr at tropopause
    94   LOGICAL, PARAMETER :: debug=.FALSE. !--- Force writefield_phy multiple outputs
     96  REAL,    PARAMETER :: gamm   =0.4   !--- Max. stretched layer sigma ratio
     97  REAL,    PARAMETER :: rho    =1.4   !--- Max tropopauses sigma ratio
     98  REAL,    PARAMETER :: o3t0   =1.E-7 !--- Nominal O3 vmr at tropopause
     99  LOGICAL, PARAMETER :: lO3Tpara=.FALSE. !--- Parametrized O3 vmr at tropopause
     100  LOGICAL, PARAMETER :: ldebug=.TRUE.!--- Force writefield_phy multiple outputs
    95101  REAL, PARAMETER :: ChemPTrMin=9.E+3 !--- Thresholds for minimum and maximum
    96102  REAL, PARAMETER :: ChemPTrMax=4.E+4 !    chemical  tropopause pressure (Pa).
    … …  
    102108!-------------------------------------------------------------------------------
    103109!
    104 SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3,             &
    105                              time_in, lon_in, lat_in, pcen_in, ptrop_ou)
     110SUBROUTINE regr_pr_time_av(fID, nam, julien, Ploc, Pre_in, Pre_ou, v3, Pgnd_ou,&
     111                                             time_in, lon_in, lat_in, Ptrp_ou)
    106112!
    107113!-------------------------------------------------------------------------------
    … …  
    118124  USE slopes_m,       ONLY: slopes
    119125  USE mod_phys_lmdz_mpi_data,       ONLY: is_mpi_root
    120   USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev
     126  USE mod_grid_phy_lmdz, ONLY: nlon=>nbp_lon, nlat=>nbp_lat, nlev_ou=>nbp_lev
    121127  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter
    122128  USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
    123129!-------------------------------------------------------------------------------
    124130! Arguments:
    125   INTEGER, INTENT(IN)  :: fID                 !--- NetCDF file ID
     131  INTEGER,           INTENT(IN) :: fID        !--- NetCDF file ID
    126132  CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
    127   REAL,    INTENT(IN)  :: julien              !--- Days since Jan 1st
    128   !--- File & LMDZ (resp. decreasing & increasing order) interfaces pressure, Pa
    129   REAL,    INTENT(IN)  :: pint_in(:)          !--- in:  file          (nlev_in+1)
    130   REAL,    INTENT(IN)  :: pint_ou(:,:)        !--- out: LMDZ     (klon,nbp_lev+1)
    131   REAL,    INTENT(OUT) :: v3(:,:,:)           !--- Regridded fld (klon,llm,n_var)
     133  REAL,              INTENT(IN) :: julien     !--- Days since Jan 1st
     134  CHARACTER(LEN=1),  INTENT(IN) :: Ploc       !--- Pressures locations
     135  !--- File/LMDZ (resp. decreasing & increasing order) pressure, Pa
     136  !    At cells centers or interfaces depending on "Ploc" keyword (C/I)
     137  REAL,    INTENT(IN)  :: Pre_in(:)           !--- in:  file      (nlev_in[+1])
     138  REAL,    INTENT(IN)  :: Pre_ou(:,:)         !--- out: LMDZ (klon,nlev_ou[+1])
     139  REAL,    INTENT(OUT) :: v3(:,:,:)           !--- Regr. fld (klon,nlev_ou,n_var)
     140  REAL,    INTENT(IN), OPTIONAL :: Pgnd_ou(:) !--- LMDZ ground pressure   (klon)
    132141  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
    133142                                              !    since Jan 1 of current year
    134143  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- File longitudes vector (klon)
    135144  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- File latitudes  vector (klon)
    136   REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- File mid-layers pres (nlev_in)
    137   REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- LMDZ tropopause pres   (klon)
     145  REAL,    INTENT(IN), OPTIONAL :: Ptrp_ou(:) !--- LMDZ tropopause pres   (klon)
    138146!-------------------------------------------------------------------------------
    139147! Local variables:
    … …  
    142150  CHARACTER(LEN=80)  :: sub
    143151  CHARACTER(LEN=320) :: str
    144   INTEGER :: vID, ncerr, n_var, ibot, ibo0, nn, itrp
    145   INTEGER :: i, nlev_in, n_dim, itop, ito0, i0
     152  INTEGER :: vID, ncerr, n_var, ibot, iout, nn
     153  INTEGER :: i, nlev_in, n_dim, itop, itrp, i0
    146154  LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
    147155  REAL    :: y_frac                           !--- Elapsed year fraction
    148156  REAL    :: alpha, beta, al                  !--- For stretching/interpolation
    149157  REAL    :: SigT_in, SigT_ou                 !--- Input and output tropopauses
    150   REAL    :: Sig_bot, Sig_top                 !--- Fully strained layer  bounds
    151   REAL    :: Sig_bo0, Sig_to0                 !--- Total strained layers bounds
    152   REAL    :: Sig_in(SIZE(pint_in))            !--- Input field sigma levels
    153   REAL    :: Sig_ou (nbp_lev+1)               !--- Output LMDZ sigma levels
    154   REAL    :: phi    (nbp_lev+1)               !--- Stretching exponent anomaly
    155   REAL    :: pstr_ou(nbp_lev+1)               !--- Stretched pressure levels
    156   REAL    :: pintou (nbp_lev+1)               !--- pint_ou in reversed order
    157   REAL    :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam))
    158   REAL    :: v2(klon,           SIZE(pint_in)-1,SIZE(nam))
    159 ! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
    160 ! v2: Field scattered to the partial "physics" horizontal grid.
     158  REAL    :: Sig_bot, Sig_top                 !--- Bounds of quasi-hat function
     159  REAL    :: Sig_bo0, Sig_to0                 !--- Lower/upper tropopauses
     160  REAL    :: Sig_in (SIZE(Pre_in))            !--- Input field sigma levels
     161  REAL    :: Sig_ou (SIZE(Pre_ou,2))          !--- Output LMDZ sigma levels
     162  REAL    :: phi    (SIZE(Pre_ou,2))          !--- Stretching exponent anomaly
     163  REAL    :: Pstr_ou(SIZE(Pre_ou,2))          !--- Stretched pressure levels
     164  REAL    :: Pres_ou(SIZE(Pre_ou,2))          !--- Pre_ou(i,:), reversed order
     165  REAL, DIMENSION(nlon, nlat) :: pg1,      pt1,      ot1
     166  REAL, DIMENSION(klon)       :: Pgnd_in,  Ptrp_in,  Otrp_in
     167  REAL, DIMENSION(klon)       :: Ptrop_ou, Pgrnd_ou
     168! * The following fields are scattered to the partial "physics" horizontal grid.
     169  REAL, POINTER :: v2(:,:,:)                  !--- Current  time ozone fields
    161170!     In the 2D file case, the values are the same for all longitudes.
    162 !     "v2(i,    k, l)" is at longitude-latitude  "xlon(i)-xlat(i)".
    163 ! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)"
    164   REAL, DIMENSION(nbp_lon, nbp_lat)   :: ps1, pt1, ot1
    165   REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
     171!     "v2(i, k, l)" is at longitude-latitude "xlon(i)-xlat(i)" and name "nam(l)"
     172! Both are:          * if Ploc=='I' in pressure interval "press_in_edg(k:k+1)"
     173!                    * if Ploc=='C' at pressure          "press_in_cen(k)"
     174  REAL, TARGET :: &
     175    v2i(klon,SIZE(Pre_in)-1,SIZE(nam)), &     !--- v2 in Ploc=='I' case
     176    v2c(klon,SIZE(Pre_in)  ,SIZE(nam))        !--- v2 in Ploc=='C' case
    166177  LOGICAL :: ll
    167178!--- For debug
    168   REAL, DIMENSION(klon)           :: ptrop_in, ptrop_ef!, dzStrain, dzStrain0
    169 ! REAL, DIMENSION(klon,nbp_lev+1) :: pstrn_ou, phii
     179  REAL, DIMENSION(klon)             :: Ptrop_in, Ptrop_ef
     180  REAL, DIMENSION(klon)             :: dzStrain, dzStrain0
     181  REAL, DIMENSION(klon,SIZE(Pre_ou,2)) :: Pstrn_ou, phii
    170182!-------------------------------------------------------------------------------
    171183  sub="regr_pr_time_av"
    172   nlev_in=SIZE(pint_in)-1
    173   CALL assert(SIZE(v3,1)==klon,          TRIM(sub)//" v3 klon")
    174   CALL assert(SIZE(v3,2)==nbp_lev,       TRIM(sub)//" v3 nbp_lev")
     184  nlev_in=SIZE(Pre_in); IF(Ploc=='I') nlev_in=nlev_in-1
     185  IF(Ploc=='I') THEN; v2 => v2i; ELSE; v2 => v2c; END IF
     186  CALL assert(SIZE(Pre_ou,1)==klon,TRIM(sub)//" Pre_ou klon")
     187  CALL assert(SIZE(v3,1)==klon,    TRIM(sub)//" v3 klon")
     188  CALL assert(SIZE(v3,2)==nlev_ou, TRIM(sub)//" v3 nlev_ou")
     189  IF(Ploc=='I') CALL assert(SIZE(Pre_ou,2)==nlev_ou+1,TRIM(sub)//" Pre_ou nlev_ou+1")
     190  IF(Ploc=='C') CALL assert(SIZE(Pre_ou,2)==nlev_ou  ,TRIM(sub)//" Pre_ou nlev_ou")
    175191  n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
    176   CALL assert(SIZE(pint_ou,1)==klon     ,TRIM(sub)//" pint_ou klon")
    177   CALL assert(SIZE(pint_ou,2)==nbp_lev+1,TRIM(sub)//" pint_ou nbp_lev+1")
    178   IF(PRESENT(lon_in))   CALL assert(SIZE(lon_in  )==klon,TRIM(sub)//" lon_in klon")
    179   IF(PRESENT(lat_in))   CALL assert(SIZE(lat_in  )==klon,TRIM(sub)//" lat_in klon")
    180   IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon")
    181   IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in )==nlev_in,TRIM(sub)//" pcen_in")
    182   lAdjTro=PRESENT(ptrop_ou)
    183   IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) &
    184   CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1)
     192  IF(PRESENT(Pgnd_ou)) CALL assert(SIZE(Pgnd_ou)==klon,TRIM(sub)//" Pgnd_ou klon")
     193  IF(PRESENT(lon_in))  CALL assert(SIZE(lon_in )==klon,TRIM(sub)//" lon_in klon")
     194  IF(PRESENT(lat_in))  CALL assert(SIZE(lat_in )==klon,TRIM(sub)//" lat_in klon")
     195  IF(PRESENT(Ptrp_ou)) CALL assert(SIZE(Ptrp_ou)==klon,TRIM(sub)//" Ptrp_ou klon")
     196  lAdjTro=PRESENT(Ptrp_ou)
     197  IF(lAdjTro) THEN
     198    IF(.NOT.PRESENT(lat_in)) &
     199      CALL abort_physic(sub, 'Missing lat_in (required if adjust_tropopause=T)', 1)
     200    IF(.NOT.PRESENT(Pgnd_ou).AND.Ploc=='C') &
     201      CALL abort_physic(sub, 'Missing ground Pr(required if adjust_tropopause=T)', 1)
     202    IF(PRESENT(Pgnd_ou)) THEN; Pgrnd_ou=Pgnd_ou; ELSE; Pgrnd_ou=Pre_ou(:,1); END IF
     203  END IF
    185204
    186205  !$OMP MASTER
    … …  
    189208    !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE
    190209    IF(lfirst) THEN
    191       lPrSurf=NF90_INQ_VARID(fID,"ps"                     ,vID)==NF90_NOERR
    192       lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
    193       lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
     210      lPrSfile=lAdjTro.AND.NF90_INQ_VARID(fID,"ps"                     ,vID)==NF90_NOERR
     211      lPrTfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
     212      lO3Tfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
    194213      CALL NF95_INQ_DIMID(fID,"time",vID)
    195214      CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
    196215      linterp=PRESENT(time_in).AND.ntim_in==14
     216      ALLOCATE(v1(nlon,nlat,nlev_in,n_var))
    197217      IF(linterp) THEN
    198         ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
    199         ALLOCATE(v1p(nbp_lon,nbp_lat,nlev_in,n_var))
    200         ALLOCATE(psm(nbp_lon,nbp_lat),psp(nbp_lon,nbp_lat))
    201         ALLOCATE(ptm(nbp_lon,nbp_lat),ptp(nbp_lon,nbp_lat))
    202         IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat))
     218        ALLOCATE(v1m(nlon,nlat,nlev_in,n_var),v1p(nlon,nlat,nlev_in,n_var))
     219        IF(lPrSfile) ALLOCATE(pgm(nlon,nlat),pgp(nlon,nlat))
     220        IF(lPrTfile) ALLOCATE(ptm(nlon,nlat),ptp(nlon,nlat))
     221        IF(lO3Tfile) ALLOCATE(otm(nlon,nlat),otp(nlon,nlat))
    203222      END IF
    204223      !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
    205       IF(lAdjTro) itrp0=locate(pcen_in,pTropUp)
    206       CALL msg(lPrSurf,'Using GROUND PRESSURE from input O3 forcing file.',sub)
    207       CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.',sub)
    208       CALL msg(lAdjTro,'o3 forcing file tropopause location uses:',sub)
    209       IF(lPrTrop) THEN
    210         CALL msg(lAdjTro,'    PRESSURE AT TROPOPAUSE from file.')
    211       ELSE IF(lO3Trop) THEN
    212         CALL msg(lAdjTro,'    O3 CONCENTRATION AT TROPOPAUSE from file.')
    213       ELSE IF(lo3tp) THEN
    214         CALL msg(lAdjTro,'    PARAMETRIZED O3 concentration at tropopause.')
    215       ELSE
    216         CALL msg(lAdjTro,'    CONSTANT O3 concentration at tropopause.')
    217       END IF
     224      IF(lAdjTro) itrp0=locate(Pre_in,pTropUp)
     225      CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.'    ,sub)
     226      CALL msg(lPrSfile,'Using GROUND PRESSURE from input O3 forcing file.',sub)
     227      CALL msg(lAdjTro ,'o3 forcing file tropopause location uses:'        ,sub)
     228      IF(lPrTfile)      THEN; str='    INPUT FILE PRESSURE'
     229      ELSE IF(lO3Tfile) THEN; str='    INPUT FILE O3 CONCENTRATION'
     230      ELSE IF(lO3Tpara) THEN; str='    PARAMETRIZED O3 concentration'
     231      ELSE;                   str='    CONSTANT O3 concentration'; END IF
     232      CALL msg(lAdjTro,TRIM(str)//' at tropopause')
    218233    END IF
    219234
    … …  
    223238    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
    224239    IF(linterp) THEN
    225       WRITE(str,'(3(a,f12.8))')'Interpolating O3 at julian day ',julien,' from '&
    226       &//'fields at times ',time_in(irec),' and ', time_in(irec+1)
     240      WRITE(str,'(a,f12.8,2(a,f5.1))')'Interpolating O3 at julian day ',julien,&
     241        ' from fields at times ',time_in(irec),' and ', time_in(irec+1)
    227242      CALL msg(.TRUE.,str,sub)
    228243      al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
    229244      v1=al*v1m+(1.-al)*v1p
    230       IF(lPrSurf) ps1=al*psm+(1.-al)*psp
    231       IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp
    232       IF(lO3Trop) ot1=al*otm+(1.-al)*otp
     245      IF(lPrSfile) pg1=al*pgm+(1.-al)*pgp
     246      IF(lPrTfile) pt1=al*ptm+(1.-al)*ptp
     247      IF(lO3Tfile) ot1=al*otm+(1.-al)*otp
    233248    END IF
    234249  END IF
    235250  !$OMP END MASTER
    236251  IF(lfirst) THEN
    237     lfirst=.FALSE.;      CALL bcast(lfirst)
    238     IF(lAdjTro)          CALL bcast(itrp0)
    239     CALL bcast(lPrTrop); CALL bcast(lPrSurf)
    240     CALL bcast(lO3Trop); CALL bcast(linterp)
     252    lfirst=.FALSE.;       CALL bcast(lfirst)
     253    IF(lAdjTro)           CALL bcast(itrp0)
     254    CALL bcast(lPrSfile); CALL bcast(lPrTfile)
     255    CALL bcast(lO3Tfile); CALL bcast(linterp)
    241256  END IF
    242257  CALL scatter2d(v1,v2)
    243   !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
    244   IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF
    245   IF(lPrTrop) CALL scatter2d(pt1,pt2)
    246   IF(lO3Trop) CALL scatter2d(ot1,ot2)
    247 
    248   !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
    249   IF(.NOT.lAdjTro) THEN
     258  IF(lPrSfile) CALL scatter2d(pg1,Pgnd_in)
     259  IF(lPrTfile) CALL scatter2d(pt1,Ptrp_in)
     260  IF(lO3Tfile) CALL scatter2d(ot1,Otrp_in)
     261  !--- No ground pressure in input file => choose it to be the one of LMDZ
     262  IF(lAdjTro.AND..NOT.lPrSfile) Pgnd_in(:)=Pgrnd_ou(:)
     263 
     264!-------------------------------------------------------------------------------
     265  IF(.NOT.lAdjTro) THEN       !--- REGRID IN PRESSURE ; NO TROPOPAUSE ADJUSTMENT
     266!-------------------------------------------------------------------------------
    250267    DO i=1,klon
    251       pintou = pint_ou(i,nbp_lev+1:1:-1)
    252       CALL regr_conserv(1,v2(i,:,:), pint_in(:), pintou(:),                    &
    253                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:), pint_in(:)))
     268      Pres_ou=Pre_ou(i,SIZE(Pre_ou,2):1:-1)   !--- pplay & paprs are decreasing
     269      IF(Ploc=='C') CALL regr_lint   (1,v2(i,:,:), LOG(Pre_in(:)),             &
     270        LOG(Pres_ou(:)), v3(i,nlev_ou:1:-1,:))
     271      IF(Ploc=='I') CALL regr_conserv(1,v2(i,:,:),     Pre_in(:) ,             &
     272            Pres_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
    254273    END DO
    255   ELSE
     274!-------------------------------------------------------------------------------
     275  ELSE                        !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT
     276!-------------------------------------------------------------------------------
    256277    y_frac=(REAL(days_elapsed)+jH_cur)/year_len
    257278
    … …  
    259280    DO i=1,klon
    260281
    261       !--- LOCAL INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS AT INTERFACES
    262       pintou=pint_ou(i,nbp_lev+1:1:-1)            !--- increasing values
    263       Sig_in(:) = [pint_in(1:nlev_in+1)/ps2(i)]   !--- increasing values
    264       Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values
     282      !--- INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS IN CURRENT COLUMN
     283      Pres_ou   = Pre_ou(i,SIZE(Pre_ou,2):1:-1)!--- pplay & paprs are decreasing
     284      Sig_in(:) = Pre_in (:)/Pgnd_in(i)            !--- increasing values
     285      Sig_ou(:) = Pres_ou(:)/Pgnd_ou(i)            !--- increasing values
    265286
    266287      !--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered
    … …  
    268289      ! ozone hole fooling the crude chemical tropopause detection algorithm.
    269290      SigT_in = get_SigTrop(i,itrp)
    270       SigT_in=MIN(SigT_in,ChemPTrMax/ps2(i))      !--- too low  value filtered
    271       SigT_in=MAX(SigT_in,ChemPTrMin/ps2(i))      !--- too high value filtered
     291      SigT_in=MIN(SigT_in,ChemPTrMax/Pgnd_in(i))   !--- too low  value filtered
     292      SigT_in=MAX(SigT_in,ChemPTrMin/Pgnd_ou(i))   !--- too high value filtered
    272293
    273294      !--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
    274       ! dynamical tropopause (especially in filaments) are conterbalanced with a
    275       ! filter ensuring it stays within a certain distance around input (file)
     295      ! dynamical tropopause (especially in filaments) are conterbalanced with
     296      ! a filter ensuring it stays within a certain distance around input (file)
    276297      ! tropopause, hence avoiding avoid a too thick stretched region ; a final
    277298      ! extra-safety filter keeps the tropopause pressure value realistic.
    278       SigT_ou = ptrop_ou(i)/ps2(i)
    279       IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho !--- too low  value w/r input
    280       IF(SigT_ou>SigT_in*rho) SigT_ou=SigT_in*rho !--- too high value w/r input
    281       SigT_ou=MIN(SigT_ou,DynPTrMax/ps2(i))       !--- too low  value filtered
    282       SigT_ou=MAX(SigT_ou,DynPTrMin/ps2(i))       !--- too high value filtered
    283       ptropou(i)=SigT_ou*ps2(i)
    284 
    285       !--- FULLY STRETCHED LAYER BOUNDS (alpha power law fully applied)
     299      SigT_ou = Ptrp_ou(i)/Pgnd_ou(i)
     300      IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho  !--- too low  value w/r input
     301      IF(SigT_ou>SigT_in*rho) SigT_ou=SigT_in*rho  !--- too high value w/r input
     302      SigT_ou=MIN(SigT_ou,DynPTrMax/Pgnd_ou(i))    !--- too low  value filtered
     303      SigT_ou=MAX(SigT_ou,DynPTrMin/Pgnd_ou(i))    !--- too high value filtered
     304      Ptrop_ou(i)=SigT_ou*Pgnd_ou(i)
     305      iout = locate(Sig_ou(:),SigT_ou)
     306
     307      !--- POWER LAW COEFFICIENT FOR TROPOPAUSES MATCHING
    286308      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
    287       Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot)
    288       Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top)
    289 
    290       !--- PARTIALLY STRETCHED LAYER BOUNDS (fraction of alpha applied)
    291       ! The used formulae ensure the first derivative stays positive.
    292       beta = LOG(Sig_top)/LOG(Sig_bot)
    293       Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
    294       Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top **  beta
    295 
    296       !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS
    297       !--- gamma<log(Sig_bo0/|alpha|) to keep Sig_bo0<1
    298       Sig_bo0 = MIN(Sig_bo0*EXP( gamm*ABS(alpha)), 0.95+(1.-0.95)*Sig_bo0)
    299       Sig_to0 =     Sig_to0*EXP(-gamm*ABS(alpha))
    300       ibo0 = locate(Sig_ou(:),Sig_bo0)
    301       ito0 = locate(Sig_ou(:),Sig_to0)
    302 
    303       !--- STRETCHING POWER LAW FUNCTION:
    304       !    0    in [0,Sig_to0] and [Sig_bo0,1] ;    1    in [Sig_bot,Sig_top]
    305       !    0->1 in [Sig_to0,Sig_top]                1->0 in [Sig_bot,Sig_bo0]
     309
     310      !--- DETERMINE STRETCHING DOMAIN UPPER AND LOWER BOUNDS
     311      Sig_bo0 = MAX(SigT_in,SigT_ou)               !--- lowest  tropopause
     312      Sig_to0 = MIN(SigT_in,SigT_ou)               !--- highest tropopause
     313      beta    = (Sig_bo0/Sig_to0)**gamm            !--- stretching exponent
     314      Sig_bot = MIN(Sig_bo0*beta,0.1*(9.+Sig_bot)) !--- must be <1
     315      ibot = locate(Sig_ou(:),Sig_bot)             !--- layer index
     316      IF(ibot-iout<2) THEN                         !--- at least one layer thick
     317        ibot=MIN(iout+2,nlev_ou); Sig_bot=Sig_ou(ibot)
     318      END IF
     319      Sig_top = Sig_to0/beta                       !--- upper bound
     320      itop = locate(Sig_ou(:),Sig_top)             !--- layer index
     321      IF(iout-itop<2) THEN                         !--- at least one layer thick
     322        itop=MAX(iout-2,1); Sig_top=Sig_ou(itop)
     323      END IF
     324
     325      !--- STRETCHING POWER LAW LOCALIZATION FUNCTION:
     326      !    0 in [0,Sig_top]    0->1 in [Sig_top,SigT_ou]
     327      !    0 in [Sig_bot,1]    1->0 in [SigT_ou, Sig_bot]
    306328      phi(:)=0.
    307       phi(itop+1:ibot) =  1.
    308       phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))&
    309                             *LOG(Sig_top)/LOG(Sig_top/Sig_to0)
    310       phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))&
    311                             *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)
     329      phi(itop+1:iout) = (1.-LOG(Sig_top)/LOG(Sig_ou(itop+1:iout)))&
     330                            *LOG(SigT_ou)/LOG(SigT_ou/Sig_top)
     331      phi(iout+1:ibot) = (1.-LOG(Sig_bot)/LOG(Sig_ou(iout+1:ibot)))&
     332                            *LOG(SigT_ou)/LOG(SigT_ou/Sig_bot)
    312333
    313334      !--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
    314       pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
     335      Pstr_ou(:) = Pres_ou(:) * Sig_ou(:)**(alpha*phi(:))
    315336
    316337      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
    317       CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:),                  &
    318                            v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:)))
     338      IF(Ploc=='C') CALL regr_lint   (1, v2(i,:,:), LOG(Pre_in(:)),            &
     339        LOG(Pstr_ou(:)), v3(i,nlev_ou:1:-1,:))
     340      IF(Ploc=='I') CALL regr_conserv(1, v2(i,:,:),     Pre_in(:) ,            &
     341            Pstr_ou(:) , v3(i,nlev_ou:1:-1,:), slopes(1,v2(i,:,:), Pre_in(:)))
    319342
    320343      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
    321       i0=nbp_lev-locate(pintou(:),ptropou(i))+1
    322       ll=check_ozone(v3(i, 1:i0     ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &
     344      i0=nlev_ou-locate(Pres_ou(:),Ptrop_ou(i))+1
     345      ll=check_ozone(v3(i, 1:i0-1   ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &
    323346                     5.E-9,3.0E-7)
    324347!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
    325       ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', &
     348      ll=check_ozone(v3(i,i0:nlev_ou,1),lon_in(i),lat_in(i),i0,'stratosphere', &
    326349                     5.E-8,1.5E-5)
    327350!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
    328351
    329       IF(debug) THEN
    330 !        dzStrain0(i)=7.*LOG(Sig_bo0/Sig_to0)
    331 !        dzStrain (i)=7.*LOG(Sig_bot/Sig_top)
    332         ptrop_in(i) = SigT_in*ps2(i)
    333 !        Pstrn_ou(i,:)= pstr_ou
    334 !        phii(i,:)=phi(:)
    335         ptrop_ef(i)=chem_tropopause(i, itrp, locate(pintou(:),PTropUp),        &
    336                            pintou(:), v3(:,nbp_lev+1:1:-1,1),o3trop=o3t0)
     352      IF(ldebug) THEN
     353        dzStrain0(i) = SIGN(7.*LOG(Sig_bo0/Sig_to0),SigT_in-SigT_ou)
     354        dzStrain (i) = SIGN(7.*LOG(Sig_bot/Sig_top),SigT_in-SigT_ou)
     355        Ptrop_in (i) = SigT_in*Pgnd_in(i)
     356        Pstrn_ou(i,:)= Pstr_ou
     357        phii(i,:)    = phi(:)
     358        Ptrop_ef(i)  = PTrop_chem(i, itrp, locate(Pres_ou(:),PTropUp),    &
     359                             Pres_ou(:), v3(:,nlev_ou:1:-1,1),o3trop=o3t0)
    337360      END IF
    338361    END DO
    339     IF(debug) THEN
    340 !      CALL writefield_phy('PreSt_ou' ,Pstrn_ou,nbp_lev+1) !--- Strained pressure
    341 !      CALL writefield_phy('dzStrain' ,dzStrain ,1)   !--- Fully & total strained
    342 !      CALL writefield_phy('dzStrain0',dzStrain0,1)   !--- regions thickness
    343 !      CALL writefield_phy('phi',phii,nbp_lev+1)      !--- Localisation function
    344       !--- Tropopauses pressures:
    345       CALL writefield_phy('PreTr_in',ptrop_in,1)     !--- Input and effective
    346       CALL writefield_phy('PreTr_ef',ptrop_ef,1)     !    chemical tropopauses
    347     END IF
    348   END IF
    349   IF(debug) THEN
    350     IF(lAdjTro) &
    351     CALL writefield_phy('PreTr_ou',ptropou,1)        !--- LMDz dyn tropopause
     362  END IF
     363  IF(ldebug.AND.lAdjTro) THEN
     364    CALL writefield_phy('PreSt_ou' ,Pstrn_ou,SIZE(Pre_ou,2)) !--- Strained Pres
     365    CALL writefield_phy('dzStrain' ,dzStrain ,1)     !--- Strained thickness
     366    CALL writefield_phy('dzStrain0',dzStrain0,1)     !--- Tropopauses distance
     367    CALL writefield_phy('phi',phii,nlev_ou)          !--- Localization function
     368    !--- Tropopauses pressures:
     369    CALL writefield_phy('PreTr_in',Ptrop_in,1)       !--- Input and effective
     370    CALL writefield_phy('PreTr_ou',Ptrop_ou,1)       !--- LMDz dyn tropopause
     371    CALL writefield_phy('PreTr_ef',Ptrop_ef,1)       !--- Effective chem tropop
     372  END IF
     373  IF(ldebug) THEN
    352374    CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
    353     CALL writefield_phy('Ozone_ou',v3(:,:,1),nbp_lev)!--- Output ozone field
    354     CALL writefield_phy('Pint_ou' ,pint_ou,1+nbp_lev)!--- LMDZ unstreched Press
     375    CALL writefield_phy('Ozone_ou',v3(:,:,1),nlev_ou)!--- Output ozone field
     376    CALL writefield_phy('Pres_ou' ,Pre_ou,SIZE(Pre_ou,2))!--- LMDZ Pressure
    355377  END IF
    356378
    … …  
    369391    CALL get_3Dfields(v1)               !--- Read ozone field(s)
    370392    IF(lAdjTro) THEN                    !--- Additional files for fields strain
    371       IF(lPrSurf) CALL get_2Dfield(ps1,"ps")
    372       IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure")
    373       IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause")
     393      IF(lPrSfile) CALL get_2Dfield(pg1,"ps")
     394      IF(lPrTfile) CALL get_2Dfield(pt1,"tropopause_air_pressure")
     395      IF(lO3Tfile) CALL get_2Dfield(ot1,"tro3_at_tropopause")
    374396    END IF
    375397  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
    … …  
    383405      CALL get_3Dfields(v1m)            !--- Read ozone field(s)
    384406      IF(lAdjTro) THEN                  !--- Additional files for fields strain
    385         IF(lPrSurf) CALL get_2Dfield(psm,"ps")
    386         IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure")
    387         IF(lO3Trop) CALL get_2Dfield(otm,"tro3_at_tropopause")
     407        IF(lPrSfile) CALL get_2Dfield(pgm,"ps")
     408        IF(lPrTfile) CALL get_2Dfield(ptm,"tropopause_air_pressure")
     409        IF(lO3Tfile) CALL get_2Dfield(otm,"tro3_at_tropopause")
    388410      END IF
    389411    ELSE                                !=== SHIFT FIELDS
    … …  
    394416      v1m=v1p                           !--- Ozone fields
    395417      IF(lAdjTro) THEN                  !--- Additional files for fields strain
    396         IF(lPrSurf) psm=psp             !--- Surface pressure
    397         IF(lPrTrop) ptm=ptp             !--- Tropopause pressure
    398         IF(lO3Trop) otm=otp             !--- Tropopause ozone
     418        IF(lPrSfile) pgm=pgp             !--- Surface pressure
     419        IF(lPrTfile) ptm=ptp             !--- Tropopause pressure
     420        IF(lO3Tfile) otm=otp             !--- Tropopause ozone
    399421      END IF
    400422    END IF
    … …  
    405427    CALL get_3Dfields(v1p)              !--- Read ozone field(s)
    406428    IF(lAdjTro) THEN                    !--- Additional files for fields strain
    407       IF(lPrSurf) CALL get_2Dfield(psp,"ps")
    408       IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure")
    409       IF(lO3Trop) CALL get_2Dfield(otp,"tro3_at_tropopause")
     429      IF(lPrSfile) CALL get_2Dfield(pgp,"ps")
     430      IF(lPrTfile) CALL get_2Dfield(ptp,"tropopause_air_pressure")
     431      IF(lO3Tfile) CALL get_2Dfield(otp,"tro3_at_tropopause")
    410432    END IF
    411433    irec=irec-1
    … …  
    438460  !--- Flip latitudes: ascending in input file, descending in "rlatu".
    439461  IF(n_dim==3) THEN
    440     v(1,:) = v(1,nbp_lat:1:-1)
    441     v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
     462    v(1,:) = v(1,nlat:1:-1)
     463    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nlon-1)  !--- Duplication
    442464  ELSE
    443     v(:,:) = v(:,nbp_lat:1:-1)
     465    v(:,:) = v(:,nlat:1:-1)
    444466  END IF
    445467
    … …  
    471493  !--- Flip latitudes: ascending in input file, descending in "rlatu".
    472494  IF(n_dim==3) THEN
    473     v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
    474     v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
     495    v(1,:,:,:) = v(1,nlat:1:-1,:,:)
     496    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1)  !--- Duplication
    475497  ELSE
    476     v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:)
     498    v(:,:,:,:) = v(:,nlat:1:-1,:,:)
    477499  END IF
    478500
    … …  
    485507!-------------------------------------------------------------------------------
    486508!
    487 FUNCTION get_SigTrop(ih,it)
    488 !
    489 !-------------------------------------------------------------------------------
    490 ! Arguments:
     509FUNCTION get_SigTrop(ih,it) RESULT(out)
     510!
     511!-------------------------------------------------------------------------------
     512! Arguments:
     513  REAL                 :: out
    491514  INTEGER, INTENT(IN)  :: ih
    492515  INTEGER, INTENT(OUT) :: it
    493   REAL                 :: get_Sigtrop
    494 !-------------------------------------------------------------------------------
    495   !--- Pressure at tropopause is read in the forcing file
    496   IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
    497     get_SigTrop=pt2(ih)/ps2(ih); RETURN
    498   END IF
    499   !--- Chemical tropopause definition is used using a particular threshold
    500   IF(lO3Trop) THEN                             !--- o3trop KNOWN FROM FILE
    501     get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,ot2(ih))
    502   ELSE IF(lo3tp) THEN                          !--- o3trop PARAMETRIZATION
    503     get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in)
    504   ELSE                                         !--- o3trop CONSTANT
    505     get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,o3t0)
    506   END IF
    507   get_SigTrop=get_SigTrop/ps2(ih)
     516!-------------------------------------------------------------------------------
     517  !--- Pressure at tropopause read from the forcing file
     518       IF(lPrTfile) THEN; out=Ptrp_in(ih)/Pgnd_in(ih); RETURN; END IF
     519
     520  !--- Chemical tropopause definition based on a particular threshold
     521       IF(lO3Tfile) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),Otrp_in(ih))
     522  ELSE IF(lO3Tpara) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1))
     523  ELSE                  ; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1),o3t0); END IF
     524  out=out/Pgnd_in(ih)
    508525
    509526END FUNCTION get_SigTrop
    … …  
    514531!-------------------------------------------------------------------------------
    515532!
    516 FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop)
     533FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out)
    517534!
    518535!-------------------------------------------------------------------------------
    … …  
    529546!-------------------------------------------------------------------------------
    530547! Arguments:
    531   REAL ::                   chem_tropopause    !--- Pressure at tropopause
     548  REAL                        :: out           !--- Pressure at tropopause
    532549  INTEGER,        INTENT(IN)  :: ih            !--- Horizontal index
    533550  INTEGER,        INTENT(OUT) :: it            !--- Index of tropopause layer
    534551  INTEGER,        INTENT(IN)  :: it0           !--- Idx: higher than tropopause
    535   REAL,           INTENT(IN)  :: pint(:)       !--- Cells-interf Pr, increasing
     552  REAL,           INTENT(IN)  :: pres(:)       !--- Pressure profile, increasing
    536553  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)
    537   REAL, OPTIONAL, INTENT(IN)  :: pcen(:)       !--- Cells-center Pr, increasing
    538554  REAL, OPTIONAL, INTENT(IN)  :: o3trop        !--- Ozone at tropopause
    539555!-------------------------------------------------------------------------------
    540556! Local variables:
    541   REAL    :: o3t                               !--- Ozone concent. at tropopause
    542   REAL    :: al                                !--- Interpolation coefficient
    543   REAL    :: coef                              !--- Coeff of latitude modulation
     557  REAL :: o3t                                  !--- Ozone concent. at tropopause
     558  REAL :: al                                   !--- Interpolation coefficient
     559  REAL :: coef                                 !--- Coeff of latitude modulation
    544560  REAL, PARAMETER :: co3(3)=[91.,28.,20.]      !--- Coeff for o3 at tropopause
    545561!-------------------------------------------------------------------------------
    … …  
    556572  it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
    557573  al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
    558   IF(PRESENT(pcen)) THEN
    559     chem_tropopause =       pcen(it)**(1.-al) * pcen(it+1)**al
    560   ELSE
    561     chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) )
    562   END IF
    563   it = locate(pint(:), chem_tropopause)        !--- pint(it)<ptrop<pint(it+1)
    564 
    565 END FUNCTION chem_tropopause
    566 !
    567 !-------------------------------------------------------------------------------
    568 
    569 
    570 !-------------------------------------------------------------------------------
    571 !
    572 FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax)
     574  IF(Ploc=='C') out =      pres(it)**(1.-al) * pres(it+1)**al
     575  IF(Ploc=='I') out = SQRT(pres(it)**(1.-al) * pres(it+2)**al *pres(it+1))
     576  it = locate(pres(:), out)                    !--- pres(it)<Ptrop<pres(it+1)
     577
     578END FUNCTION PTrop_chem
     579!
     580!-------------------------------------------------------------------------------
     581
     582
     583!-------------------------------------------------------------------------------
     584!
     585FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out)
    573586!
    574587!-------------------------------------------------------------------------------
    … …  
    576589!-------------------------------------------------------------------------------
    577590! Arguments:
    578   LOGICAL                      :: check_ozone      !--- .T. => some wrong values
     591  LOGICAL                      :: out          !--- .T. => some wrong values
    579592  REAL,             INTENT(IN) :: o3col(:), lon, lat
    580593  INTEGER,          INTENT(IN) :: ilev0
    … …  
    592605  lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
    593606  lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
    594   check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN
    595   IF(prt_level<1) RETURN
     607  out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN
    596608
    597609  !--- SOME TOO LOW VALUES FOUND
Note: See TracChangeset for help on using the changeset viewer.