Ignore:
Timestamp:
Nov 5, 2018, 3:24:59 PM (6 years ago)
Author:
Laurent Fairhead
Message:

Undoing merge with trunk (r3356) to properly register Yann's latest modifications

Location:
LMDZ6/branches/DYNAMICO-conv
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/DYNAMICO-conv

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_pr_time_av_m.F90

    r3356 r3411  
    55  USE mod_phys_lmdz_transfert_para, ONLY: bcast
    66  USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank
    7   USE print_control_mod,  ONLY: prt_level
    87  IMPLICIT NONE
    98
    … …  
    1615!    to the LMDZ vertical grid.
    1716!  * the forcing fields are stretched if the following arguments are present:
    18 !     - "lat_in":  input file latitudes.
    19 !     - "Ptrp_ou": target grid (LMDZ) tropopause pressure.
     17!     - "lat_in":   input file latitudes.
     18!     - "pcen_in":  input file cells center pressure levels.
     19!     - "ptrop_ou": target grid (LMDZ) tropopause pressure.
    2020!   so that the tropopause is shifted to the position of the one of LMDZ.
    2121!  Note that the ozone quantity conservation is not ensured.
    … …  
    3434!    except that the latitudes are in ascending order in the input file.
    3535!  * We assume that all the input fields have the same coordinates.
    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.
     36!  * The target vertical LMDZ grid is the grid of layer boundaries.
     37!  * Regridding in pressure is conservative, second order.
    4338!  * All the fields are regridded as a single multi-dimensional array, so it
    4439!    saves CPU time to call this procedure once for several NetCDF variables
    … …  
    4742!    1) read from the file if "tropopause_air_pressure" field is available.
    4843!    2) computed using "tro3" and "tro3_at_tropopause' (if available).
    49 !    3) computed using "tro3" and a fixed threshold otherwise, constant or
    50 !    determined using an empirical three parameters law:
     44!    3) computed using "tro3" and a fixed threshold otherwise, determined using
     45!    an empirical three parameters law:
    5146!         o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3)
    5247!       => co1 and co2 are in ppbV, and co3 in degrees.
    … …  
    5651!  * Fields with suffix "m"/"p" are at the closest records earlier/later than
    5752!  the mid-julian day "julien", on the global "dynamics" horizontal grid.
    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)]".
     53!  * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
     54!    interval "pint_in(k:k+1)]" and variable "nam(l)".
    6055!  * In the 2D file case, the values are the same for all longitudes.
    6156!  * The tropopause correction works like this: the input fields (file) are
    6257!  interpolated on output (LMDZ) pressure field, which is streched using a power
    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)
     58!  law in a limited zone made of 3 layers:
     59!    1) between the two tropopauses (file and LMDZ)
     60!    2) in an upper and a lower transitions layers
    6661!  The stretching function has the following form:
    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.
     62!    Sigma_str = Sigma^(1+alpha), with alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
     63!  This value shifts the file tropopause to the height of the one of LMDZ.
     64!  The stretching is fully applied in the central zone only, and only partially
     65!  in the transitions zones, thick enough to guarantee a growing stretched
     66!  pressure field. The ponderation function for alpha to modulate the stretching
     67!  is constant equal to 1 in the central layer, and quasi-linear (from 1 to 0)
     68!  in the transition layers (from 1 to 0 ; in fact: sections of 1/log function),
     69!  making this localisation function quasi-trapezoidal.
    7570!
    7671! * The following fields are on the global "dynamics" grid, as read from files:
    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
     72  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
     73  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next ozone fields
     74  REAL,    SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure
    8275  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
    8376  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
    … …  
    8881!      * for monthly input files: julien is in [time_in(irec),time_in(irec+1)]
    8982  LOGICAL, SAVE :: linterp                         !--- Interpolation in time
    90   LOGICAL, SAVE :: lPrSfile                        !--- Surface pressure flag
    91   LOGICAL, SAVE :: lPrTfile                        !--- Tropopause pressure flag
    92   LOGICAL, SAVE :: lO3Tfile                        !--- Tropopause ozone flag
     83  LOGICAL, SAVE :: lPrSurf                         !--- Surface pressure flag
     84  LOGICAL, SAVE :: lPrTrop                         !--- Tropopause pressure flag
     85  LOGICAL, SAVE :: lO3Trop                         !--- Tropopause ozone flag
    9386  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
    9487!$OMP THREADPRIVATE(lfirst)
    95   REAL,    PARAMETER :: pTropUp=9.E+3 !--- Value  <  tropopause pressure (Pa)
    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=.FALSE.!--- Force writefield_phy multiple outputs
    101   REAL, PARAMETER :: ChemPTrMin=9.E+3 !--- Thresholds for minimum and maximum
    102   REAL, PARAMETER :: ChemPTrMax=4.E+4 !    chemical  tropopause pressure (Pa).
    103   REAL, PARAMETER :: DynPTrMin =8.E+3 !--- Thresholds for minimum and maximum
    104   REAL, PARAMETER :: DynPTrMax =4.E+4 !    dynamical tropopause pressure (Pa).
     88  REAL,    PARAMETER :: pTropUp=5.E+3 !--- Value < tropopause pressure (Pa)
     89  REAL,    PARAMETER :: gamm = 0.4    !--- Relative thickness of transitions
     90  REAL,    PARAMETER :: rho  = 1.3    !--- Max tropopauses sigma ratio
     91  REAL,    PARAMETER :: o3t0 = 1.E-7  !--- Nominal O3 vmr at tropopause
     92  LOGICAL, PARAMETER :: lo3tp=.FALSE. !--- Use parametrized O3 vmr at tropopause
    10593
    10694CONTAINS
    … …  
    10896!-------------------------------------------------------------------------------
    10997!
    110 SUBROUTINE regr_pr_time_av(fID, nam, julien, Ploc, Pre_in, Pre_ou, v3, Pgnd_ou,&
    111                                              time_in, lon_in, lat_in, Ptrp_ou)
     98SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3,             &
     99                             time_in, lon_in, lat_in, pcen_in, ptrop_ou)
    112100!
    113101!-------------------------------------------------------------------------------
    … …  
    124112  USE slopes_m,       ONLY: slopes
    125113  USE mod_phys_lmdz_mpi_data,       ONLY: is_mpi_root
    126   USE mod_grid_phy_lmdz, ONLY: nlon=>nbp_lon, nlat=>nbp_lat, nlev_ou=>nbp_lev
     114  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev
    127115  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter
    128116  USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
    129117!-------------------------------------------------------------------------------
    130118! Arguments:
    131   INTEGER,           INTENT(IN) :: fID        !--- NetCDF file ID
     119  INTEGER, INTENT(IN)  :: fID           !--- NetCDF file ID
    132120  CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
    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)
     121  REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
     122  REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces file Pres, Pa, ascending
     123  REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ pressure, Pa (g,nbp_lev+1)
     124  REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
    141125  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
    142126                                              !    since Jan 1 of current year
    143   REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- File longitudes vector (klon)
    144   REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- File latitudes  vector (klon)
    145   REAL,    INTENT(IN), OPTIONAL :: Ptrp_ou(:) !--- LMDZ tropopause pres   (klon)
     127  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- Input/output longitudes vector
     128  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
     129  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure
     130  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pres (klon)
    146131!-------------------------------------------------------------------------------
    147132! Local variables:
    … …  
    150135  CHARACTER(LEN=80)  :: sub
    151136  CHARACTER(LEN=320) :: str
    152   INTEGER :: vID, ncerr, n_var, ibot, iout, nn
    153   INTEGER :: i, nlev_in, n_dim, itop, itrp, i0
     137  INTEGER :: vID, ncerr, n_var, ibot, ibo0, nn, itrp
     138  INTEGER :: i, nlev_in, n_dim, itop, ito0, i0
    154139  LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
    155140  REAL    :: y_frac                           !--- Elapsed year fraction
    156141  REAL    :: alpha, beta, al                  !--- For stretching/interpolation
    157142  REAL    :: SigT_in, SigT_ou                 !--- Input and output tropopauses
    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
     143  REAL    :: Sig_bot, Sig_top                 !--- Fully strained layer  bounds
     144  REAL    :: Sig_bo0, Sig_to0                 !--- Total strained layers bounds
     145  REAL    :: Sig_in(SIZE(pint_in))            !--- Input field sigma levels
     146  REAL    :: Sig_ou (nbp_lev+1)               !--- Output LMDZ sigma levels
     147  REAL    :: phi    (nbp_lev+1)               !--- Stretching exponent anomaly
     148  REAL    :: pstr_ou(nbp_lev+1)               !--- Stretched pressure levels
     149  REAL    :: pintou (nbp_lev+1)               !--- pint_ou in reversed order
     150  REAL    :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam))
     151  REAL    :: v2(klon,           SIZE(pint_in)-1,SIZE(nam))
     152! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.
     153! v2: Field scattered to the partial "physics" horizontal grid.
    170154!     In the 2D file case, the values are the same for all longitudes.
    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
     155!     "v2(i,    k, l)" is at longitude-latitude  "xlon(i)-xlat(i)".
     156! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)"
     157  REAL, DIMENSION(nbp_lon, nbp_lat)   :: ps1, pt1, ot1
     158  REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
    177159  LOGICAL :: ll
    178 !--- For debug
    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
    182160!-------------------------------------------------------------------------------
    183161  sub="regr_pr_time_av"
    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")
     162  nlev_in=SIZE(pint_in)-1
     163  CALL assert(SIZE(v3,1)==klon,          TRIM(sub)//" v3 klon")
     164  CALL assert(SIZE(v3,2)==nbp_lev,       TRIM(sub)//" v3 nbp_lev")
    191165  n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
    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
     166  CALL assert(SIZE(pint_ou,1)==klon     ,TRIM(sub)//" pint_ou klon")
     167  CALL assert(SIZE(pint_ou,2)==nbp_lev+1,TRIM(sub)//" pint_ou nbp_lev+1")
     168  IF(PRESENT(lon_in))   CALL assert(SIZE(lon_in  )==klon,TRIM(sub)//" lon_in klon")
     169  IF(PRESENT(lat_in))   CALL assert(SIZE(lat_in  )==klon,TRIM(sub)//" lat_in klon")
     170  IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon")
     171  IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in )==nlev_in,TRIM(sub)//" pcen_in")
     172  lAdjTro=PRESENT(ptrop_ou)
     173  IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) &
     174  CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1)
    204175
    205176  !$OMP MASTER
    … …  
    208179    !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE
    209180    IF(lfirst) THEN
    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
     181      lPrSurf=NF90_INQ_VARID(fID,"ps"                     ,vID)==NF90_NOERR
     182      lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
     183      lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
    213184      CALL NF95_INQ_DIMID(fID,"time",vID)
    214185      CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
    215186      linterp=PRESENT(time_in).AND.ntim_in==14
    216       ALLOCATE(v1(nlon,nlat,nlev_in,n_var))
    217187      IF(linterp) THEN
    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))
     188        ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
     189        ALLOCATE(v1p(nbp_lon,nbp_lat,nlev_in,n_var))
     190        ALLOCATE(psm(nbp_lon,nbp_lat),psp(nbp_lon,nbp_lat))
     191        ALLOCATE(ptm(nbp_lon,nbp_lat),ptp(nbp_lon,nbp_lat))
     192        IF(lO3Trop) ALLOCATE(otm(nbp_lon,nbp_lat),otp(nbp_lon,nbp_lat))
    222193      END IF
    223194      !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
    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')
     195      IF(lAdjTro) itrp0=locate(pcen_in,pTropUp)
     196      CALL msg(lPrSurf,'Using GROUND PRESSURE from input O3 forcing file.',sub)
     197      CALL msg(linterp,'Monthly O3 files => ONLINE TIME INTERPOLATION.',sub)
     198      CALL msg(lAdjTro,'o3 forcing file tropopause location uses:',sub)
     199      IF(lPrTrop) THEN
     200        CALL msg(lAdjTro,'    PRESSURE AT TROPOPAUSE from file.')
     201      ELSE IF(lO3Trop) THEN
     202        CALL msg(lAdjTro,'    O3 CONCENTRATION AT TROPOPAUSE from file.')
     203      ELSE IF(lo3tp) THEN
     204        CALL msg(lAdjTro,'    PARAMETRIZED O3 concentration at tropopause.')
     205      ELSE
     206        CALL msg(lAdjTro,'    CONSTANT O3 concentration at tropopause.')
     207      END IF
    233208    END IF
    234209
    … …  
    238213    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
    239214    IF(linterp) THEN
    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)
     215      WRITE(str,'(3(a,f12.8))')'Interpolating O3 at julian day ',julien,' from '&
     216      &//'fields at times ',time_in(irec),' and ', time_in(irec+1)
    242217      CALL msg(.TRUE.,str,sub)
    243218      al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
    244219      v1=al*v1m+(1.-al)*v1p
    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
     220      IF(lPrSurf) ps1=al*psm+(1.-al)*psp
     221      IF(lPrTrop) pt1=al*ptm+(1.-al)*ptp
     222      IF(lO3Trop) ot1=al*otm+(1.-al)*otp
    248223    END IF
    249224  END IF
    250225  !$OMP END MASTER
    251226  IF(lfirst) THEN
    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)
     227    lfirst=.FALSE.;      CALL bcast(lfirst)
     228    IF(lAdjTro)          CALL bcast(itrp0)
     229    CALL bcast(lPrTrop); CALL bcast(lPrSurf)
     230    CALL bcast(lO3Trop); CALL bcast(linterp)
    256231  END IF
    257232  CALL scatter2d(v1,v2)
    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 !-------------------------------------------------------------------------------
     233  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
     234  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF
     235  IF(lPrTrop) CALL scatter2d(pt1,pt2)
     236  IF(lO3Trop) CALL scatter2d(ot1,ot2)
     237
     238  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
     239  IF(.NOT.lAdjTro) THEN
    267240    DO i=1,klon
    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(:)))
     241      pintou = pint_ou(i,nbp_lev+1:1:-1)
     242      CALL regr_conserv(1,v2(i,:,:), pint_in(:), pintou(:),                    &
     243                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:), pint_in(:)))
    273244    END DO
    274 !-------------------------------------------------------------------------------
    275   ELSE                        !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT
    276 !-------------------------------------------------------------------------------
     245  ELSE
    277246    y_frac=(REAL(days_elapsed)+jH_cur)/year_len
    278247
    … …  
    280249    DO i=1,klon
    281250
    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
    286 
    287       !--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered
    288       ! to keep tropopause pressure realistic ; high values are usually due to
    289       ! ozone hole fooling the crude chemical tropopause detection algorithm.
     251      !--- LOCAL INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS AT INTERFACES
     252      pintou=pint_ou(i,nbp_lev+1:1:-1)            !--- increasing values
     253      Sig_in(:) = [pint_in(1:nlev_in+1)/ps2(i)]   !--- increasing values
     254      Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values
     255
     256      !--- INPUT (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE
    290257      SigT_in = get_SigTrop(i,itrp)
    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
    293 
    294       !--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
    295       ! dynamical tropopause (especially in filaments) are conterbalanced with
    296       ! a filter ensuring it stays within a certain distance around input (file)
    297       ! tropopause, hence avoiding avoid a too thick stretched region ; a final
    298       ! extra-safety filter keeps the tropopause pressure value realistic.
    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
     258      SigT_ou = ptrop_ou(i)/ps2(i)
     259
     260      !--- AVOID THE FILAMENTS WHICH WOULD NEED A VERY THICK STRETCHED REGION
     261      IF(SigT_ou>SigT_in*rho) SigT_ou = SigT_in*rho
     262      IF(SigT_ou<SigT_in/rho) SigT_ou = SigT_in/rho
     263      ptropou(i)=SigT_ou*ps2(i)
     264
     265      !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW
    308266      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
    309267
    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_bo0)) !--- 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]
     268      !--- FULLY STRETCHED LAYER BOUNDS (FILE AND MODEL TROPOPAUSES)
     269      Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot)
     270      Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top)
     271
     272      !--- PARTIALLY STRETCHED LAYER BOUNDS, ENSURING >0 DERIVATIVE
     273      beta = LOG(Sig_top)/LOG(Sig_bot)
     274      Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
     275      Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top **  beta
     276
     277      !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS
     278      !--- gamma<log(Sig_bo0/|alpha|) to keep Sig_bo0<1
     279      Sig_bo0 = MIN(Sig_bo0*EXP( gamm*ABS(alpha)), 0.95+(1.-0.95)*Sig_bo0)
     280      Sig_to0 =     Sig_to0*EXP(-gamm*ABS(alpha))
     281      ibo0 = locate(Sig_ou(:),Sig_bo0)
     282      ito0 = locate(Sig_ou(:),Sig_to0)
     283
     284      !--- FUNCTION FOR STRETCHING LOCALISATION
     285      !    0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1
    328286      phi(:)=0.
    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)
    333 
    334       !--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
    335       Pstr_ou(:) = Pres_ou(:) * Sig_ou(:)**(alpha*phi(:))
     287      phi(itop+1:ibot) =  1.
     288      phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))&
     289                            *LOG(Sig_top)/LOG(Sig_top/Sig_to0)
     290      phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))&
     291                            *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)
     292
     293      !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
     294      pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
    336295
    337296      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
    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(:)))
     297      CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:),                  &
     298                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:)))
    342299
    343300      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
    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',  &
     301      i0=nbp_lev-locate(pintou(:),ptropou(i))+1
     302      ll=check_ozone(v3(i, 1:i0     ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &
    346303                     5.E-9,3.0E-7)
    347304!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
    348       ll=check_ozone(v3(i,i0:nlev_ou,1),lon_in(i),lat_in(i),i0,'stratosphere', &
     305      ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', &
    349306                     5.E-8,1.5E-5)
    350307!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
    351308
    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)
    360       END IF
    361309    END DO
    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
    374     CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
    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
    377310  END IF
    378311
    … …  
    391324    CALL get_3Dfields(v1)               !--- Read ozone field(s)
    392325    IF(lAdjTro) THEN                    !--- Additional files for fields strain
    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")
     326      IF(lPrSurf) CALL get_2Dfield(ps1,"ps")
     327      IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure")
     328      IF(lO3Trop) CALL get_2Dfield(ot1,"tro3_at_tropopause")
    396329    END IF
    397330  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
    … …  
    405338      CALL get_3Dfields(v1m)            !--- Read ozone field(s)
    406339      IF(lAdjTro) THEN                  !--- Additional files for fields strain
    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")
     340        IF(lPrSurf) CALL get_2Dfield(psm,"ps")
     341        IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure")
     342        IF(lO3Trop) CALL get_2Dfield(otm,"tro3_at_tropopause")
    410343      END IF
    411344    ELSE                                !=== SHIFT FIELDS
    … …  
    416349      v1m=v1p                           !--- Ozone fields
    417350      IF(lAdjTro) THEN                  !--- Additional files for fields strain
    418         IF(lPrSfile) pgm=pgp             !--- Surface pressure
    419         IF(lPrTfile) ptm=ptp             !--- Tropopause pressure
    420         IF(lO3Tfile) otm=otp             !--- Tropopause ozone
     351        IF(lPrSurf) psm=psp             !--- Surface pressure
     352        IF(lPrTrop) ptm=ptp             !--- Tropopause pressure
     353        IF(lO3Trop) otm=otp             !--- Tropopause ozone
    421354      END IF
    422355    END IF
    … …  
    427360    CALL get_3Dfields(v1p)              !--- Read ozone field(s)
    428361    IF(lAdjTro) THEN                    !--- Additional files for fields strain
    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")
     362      IF(lPrSurf) CALL get_2Dfield(psp,"ps")
     363      IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure")
     364      IF(lO3Trop) CALL get_2Dfield(otp,"tro3_at_tropopause")
    432365    END IF
    433366    irec=irec-1
    … …  
    460393  !--- Flip latitudes: ascending in input file, descending in "rlatu".
    461394  IF(n_dim==3) THEN
    462     v(1,:) = v(1,nlat:1:-1)
    463     v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nlon-1)  !--- Duplication
     395    v(1,:) = v(1,nbp_lat:1:-1)
     396    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
    464397  ELSE
    465     v(:,:) = v(:,nlat:1:-1)
     398    v(:,:) = v(:,nbp_lat:1:-1)
    466399  END IF
    467400
    … …  
    493426  !--- Flip latitudes: ascending in input file, descending in "rlatu".
    494427  IF(n_dim==3) THEN
    495     v(1,:,:,:) = v(1,nlat:1:-1,:,:)
    496     v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nlon-1)  !--- Duplication
     428    v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
     429    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
    497430  ELSE
    498     v(:,:,:,:) = v(:,nlat:1:-1,:,:)
     431    v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:)
    499432  END IF
    500433
    … …  
    507440!-------------------------------------------------------------------------------
    508441!
    509 FUNCTION get_SigTrop(ih,it) RESULT(out)
    510 !
    511 !-------------------------------------------------------------------------------
    512 ! Arguments:
    513   REAL                 :: out
     442FUNCTION get_SigTrop(ih,it)
     443!
     444!-------------------------------------------------------------------------------
     445! Arguments:
    514446  INTEGER, INTENT(IN)  :: ih
    515447  INTEGER, INTENT(OUT) :: it
    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)
     448  REAL                 :: get_Sigtrop
     449!-------------------------------------------------------------------------------
     450  !--- Pressure at tropopause is read in the forcing file
     451  IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
     452    get_SigTrop=pt2(ih)/ps2(ih); RETURN
     453  END IF
     454  !--- Chemical tropopause definition is used using a particular threshold
     455  IF(lO3Trop) THEN                             !--- o3trop KNOWN FROM FILE
     456    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,ot2(ih))
     457  ELSE IF(lo3tp) THEN                          !--- o3trop PARAMETRIZATION
     458    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in)
     459  ELSE                                         !--- o3trop CONSTANT
     460    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,o3t0)
     461  END IF
     462  get_SigTrop=get_SigTrop/ps2(ih)
    525463
    526464END FUNCTION get_SigTrop
    … …  
    531469!-------------------------------------------------------------------------------
    532470!
    533 FUNCTION PTrop_chem(ih,it,it0,pres,o3,o3trop) RESULT(out)
     471FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop)
    534472!
    535473!-------------------------------------------------------------------------------
    … …  
    546484!-------------------------------------------------------------------------------
    547485! Arguments:
    548   REAL                        :: out           !--- Pressure at tropopause
     486  REAL ::                   chem_tropopause    !--- Pressure at tropopause
    549487  INTEGER,        INTENT(IN)  :: ih            !--- Horizontal index
    550488  INTEGER,        INTENT(OUT) :: it            !--- Index of tropopause layer
    551489  INTEGER,        INTENT(IN)  :: it0           !--- Idx: higher than tropopause
    552   REAL,           INTENT(IN)  :: pres(:)       !--- Pressure profile, increasing
     490  REAL,           INTENT(IN)  :: pint(:)       !--- Cells-interf Pr, increasing
    553491  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)
     492  REAL, OPTIONAL, INTENT(IN)  :: pcen(:)       !--- Cells-center Pr, increasing
    554493  REAL, OPTIONAL, INTENT(IN)  :: o3trop        !--- Ozone at tropopause
    555494!-------------------------------------------------------------------------------
    556495! Local variables:
    557   REAL :: o3t                                  !--- Ozone concent. at tropopause
    558   REAL :: al                                   !--- Interpolation coefficient
    559   REAL :: coef                                 !--- Coeff of latitude modulation
     496  REAL    :: o3t                               !--- Ozone concent. at tropopause
     497  REAL    :: al                                !--- Interpolation coefficient
     498  REAL    :: coef                              !--- Coeff of latitude modulation
    560499  REAL, PARAMETER :: co3(3)=[91.,28.,20.]      !--- Coeff for o3 at tropopause
    561500!-------------------------------------------------------------------------------
    … …  
    572511  it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
    573512  al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
    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 
    578 END FUNCTION PTrop_chem
    579 !
    580 !-------------------------------------------------------------------------------
    581 
    582 
    583 !-------------------------------------------------------------------------------
    584 !
    585 FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax) RESULT(out)
     513  IF(PRESENT(pcen)) THEN
     514    chem_tropopause =       pcen(it)**(1.-al) * pcen(it+1)**al
     515  ELSE
     516    chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) )
     517  END IF
     518  it = locate(pint(:), chem_tropopause)        !--- pint(it)<ptrop<pint(it+1)
     519
     520END FUNCTION chem_tropopause
     521!
     522!-------------------------------------------------------------------------------
     523
     524
     525!-------------------------------------------------------------------------------
     526!
     527FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax)
    586528!
    587529!-------------------------------------------------------------------------------
    … …  
    589531!-------------------------------------------------------------------------------
    590532! Arguments:
    591   LOGICAL                      :: out          !--- .T. => some wrong values
     533  LOGICAL                      :: check_ozone      !--- .T. => some wrong values
    592534  REAL,             INTENT(IN) :: o3col(:), lon, lat
    593535  INTEGER,          INTENT(IN) :: ilev0
    … …  
    605547  lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
    606548  lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
    607   out=lmin.OR.lmax; IF(.NOT.out.OR.prt_level>100) RETURN
     549  check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN
    608550
    609551  !--- SOME TOO LOW VALUES FOUND
Note: See TracChangeset for help on using the changeset viewer.