source: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_pr_time_av_m.F90 @ 5017

! $Id$
MODULE regr_pr_time_av_m

  USE write_field_phy
  USE mod_phys_lmdz_transfert_para, ONLY: bcast
  USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank
  IMPLICIT NONE

!-------------------------------------------------------------------------------
! Purpose: Regrid pressure with an averaging method. Operations done:
!  * on the root process: read a NetCDF 3D or 4D field at a given day.
!  * pack the fields to the LMDZ "horizontal "physics" grid and scatter
!    to all threads of all processes;
!  * in all the threads of all the processes, regrid the fields in pressure
!    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.
!   so that the tropopause is shifted to the position of the one of LMDZ.
!  Note that the ozone quantity conservation is not ensured.
!-------------------------------------------------------------------------------
! Initial routine: regr_pr_av_m module (L. Guez).
! Present version: David Cugnet ; corresponding additions:
!    - time interpolation
!    - 3D compliant
!    - vertical stretching of the field to allow tropopause matching between
!    input field and current lmdz field.
!-------------------------------------------------------------------------------
! Remarks:
!  * 3D fields are zonal means, with dimensions (latitude, pressure, julian day)
!  * 4D fields have the dimensions:  (longitude, latitude, pressure, julian day)
!  * All the fields are already on the horizontal grid (rlatu) or (rlatu,rlonv),
!    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.
!  * 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
!    rather than several times, each time for a single NetCDF variable.
!  * The input file pressure at tropopause can be (in decreasing priority):
!    1) read from the file if "tropopause_air_pressure" field is available.
!    2) computed using "tro3" and "tro3_at_tropopause' (if available).
!    3) computed using "tro3" and a fixed threshold otherwise, determined using
!    an empirical three parameters law:
!         o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3)
!       => co1 and co2 are in ppbV, and co3 in degrees.
!       => co3 allow a smooth transition between north and south hemispheres.
!  * If available, the field "ps" (input file pressure at surface) is used.
!    If not, the current LMDZ ground pressure is taken instead.
!  * Fields with suffix "m"/"p" are at the closest records earlier/later than
!  the mid-julian day "julien", on the global "dynamics" horizontal grid.
!  * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
!    interval "pint_in(k:k+1)]" and variable "nam(l)".
!  * In the 2D file case, the values are the same for all longitudes.
!  * The tropopause correction works like this: the input fields (file) are
!  interpolated on output (LMDZ) pressure field, which is streched using a power
!  law in a limited zone made of 3 layers:
!    1) between the two tropopauses (file and LMDZ)
!    2) in an upper and a lower transitions layers
!  The stretching function has the following form:
!    Sigma_str = Sigma^(1+alpha), with alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
!  This value shifts the file tropopause to the height of the one of LMDZ.
!  The stretching is fully applied in the central zone only, and only partially
!  in the transitions zones, thick enough to guarantee a growing stretched
!  pressure field. The ponderation function for alpha to modulate the stretching
!  is constant equal to 1 in the central layer, and quasi-linear (from 1 to 0)
!  in the transition layers (from 1 to 0 ; in fact: sections of 1/log function),
!  making this localisation function quasi-trapezoidal.
!
! * The following fields are on the global "dynamics" grid, as read from files:
  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next ozone fields
  REAL,    SAVE, ALLOCATABLE :: psm(:,:), psp(:,:) !--- Surface pressure
  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
  INTEGER, SAVE :: ntim_in                         !--- Records nb in input file
  INTEGER, SAVE :: itrp0                           !--- idx above chem tropop.
  INTEGER, SAVE :: irec                            !--- Current time index
!      * for daily   input files: current julian day number
!      * 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 :: lfirst=.TRUE.                   !--- First call flag
!$OMP THREADPRIVATE(lfirst)
  REAL,    PARAMETER :: pTropUp=5.E+3 !--- Value < tropopause pressure (Pa)
  REAL,    PARAMETER :: gamm = 0.4    !--- Relative thickness of transitions
  REAL,    PARAMETER :: rho  = 1.3    !--- Max tropopauses sigma ratio
  REAL,    PARAMETER :: o3t0 = 1.E-7  !--- Nominal O3 vmr at tropopause
  LOGICAL, PARAMETER :: lo3tp=.FALSE. !--- Use parametrized O3 vmr at tropopause

CONTAINS

!-------------------------------------------------------------------------------
!
SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3,             &
                             time_in, lon_in, lat_in, pcen_in, ptrop_ou)
!
!-------------------------------------------------------------------------------
  USE dimphy,         ONLY: klon
  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, &
                            NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION
  USE netcdf,         ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR
  USE assert_m,       ONLY: assert
  USE assert_eq_m,    ONLY: assert_eq
!!  USE comvert_mod,    ONLY: scaleheight
  USE interpolation,  ONLY: locate
  USE regr_conserv_m, ONLY: regr_conserv
  USE regr_lint_m,    ONLY: regr_lint
  USE slopes_m,       ONLY: slopes
  USE mod_phys_lmdz_para,           ONLY: is_mpi_root,is_master
  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid_type, unstructured
  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter, gather
  USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
  USE geometry_mod,                 ONLY: ind_cell_glo
!-------------------------------------------------------------------------------
! Arguments:
  INTEGER, INTENT(IN)  :: fID           !--- NetCDF file ID
  CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
  REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
  REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces file Pres, Pa, ascending
  REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ pressure, Pa (g,nbp_lev+1)
  REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
                                              !    since Jan 1 of current year
  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- Input/output longitudes vector
  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure
  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pres (klon)
!-------------------------------------------------------------------------------
! Local variables:
  include "clesphys.h"
  include "YOMCST.h"
  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
  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.
!     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
  INTEGER,ALLOCATABLE :: ind_cell_glo_glo(:)
  LOGICAL :: ll
!-------------------------------------------------------------------------------
  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")
  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)

  !$OMP MASTER
  IF(is_mpi_root) THEN

    !=== 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
      CALL NF95_INQ_DIMID(fID,"time",vID)
      CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
      linterp=PRESENT(time_in).AND.ntim_in==14
      IF(linterp) THEN
        ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
        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))
      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
    END IF

    !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES)
    CALL update_fields()

    !=== 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)
      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
    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)
  END IF
  
  IF (is_master) THEN 
    ALLOCATE(ind_cell_glo_glo(klon_glo))
  ELSE
    ALLOCATE(ind_cell_glo_glo(0))
  ENDIF
  CALL gather(ind_cell_glo,ind_cell_glo_glo)
  IF (is_master .AND. grid_type==unstructured) v1(:,:,:,:)=v1(:,ind_cell_glo_glo(:),:,:)
  
  CALL scatter2d(v1,v2)

  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
  IF(lPrSurf) THEN
    IF (is_master .AND. grid_type==unstructured) ps1(:,:)=ps1(:,ind_cell_glo_glo(:))
    CALL scatter2d(ps1,ps2)
  ELSE
   ps2=pint_ou(:,1)
  END IF

  IF(lPrTrop) THEN
    IF (is_master .AND. grid_type==unstructured) pt1(:,:)=pt1(:,ind_cell_glo_glo(:))
    CALL scatter2d(pt1,pt2)
  ENDIF

  IF(lO3Trop) THEN
    IF (is_master .AND. grid_type==unstructured) ot1(:,:)=ot1(:,ind_cell_glo_glo(:))
    CALL scatter2d(ot1,ot2)
  ENDIF
  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
  IF(.NOT.lAdjTro) THEN
    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(:)))
    END DO
  ELSE
    y_frac=(REAL(days_elapsed)+jH_cur)/year_len

    !--- OUTPUT SIGMA LEVELS
    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 (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE
      SigT_in = get_SigTrop(i,itrp)
      SigT_ou = ptrop_ou(i)/ps2(i)

      !--- AVOID THE FILAMENTS WHICH WOULD NEED A VERY THICK STRETCHED REGION
      IF(SigT_ou>SigT_in*rho) SigT_ou = SigT_in*rho
      IF(SigT_ou<SigT_in/rho) SigT_ou = SigT_in/rho
      ptropou(i)=SigT_ou*ps2(i)

      !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW
      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)

      !--- FULLY STRETCHED LAYER BOUNDS (FILE AND MODEL TROPOPAUSES)
      Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot)
      Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top)

      !--- PARTIALLY STRETCHED LAYER BOUNDS, ENSURING >0 DERIVATIVE
      beta = LOG(Sig_top)/LOG(Sig_bot)
      Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
      Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top **  beta

      !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS
      !--- gamma<log(Sig_bo0/|alpha|) to keep Sig_bo0<1
      Sig_bo0 = MIN(Sig_bo0*EXP( gamm*ABS(alpha)), 0.95+(1.-0.95)*Sig_bo0)
      Sig_to0 =     Sig_to0*EXP(-gamm*ABS(alpha))
      ibo0 = locate(Sig_ou(:),Sig_bo0)
      ito0 = locate(Sig_ou(:),Sig_to0)

      !--- FUNCTION FOR STRETCHING LOCALISATION
      !    0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1
      phi(:)=0.
      phi(itop+1:ibot) =  1.
      phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))&
                            *LOG(Sig_top)/LOG(Sig_top/Sig_to0)
      phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))&
                            *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)

      !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
      pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))

      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
      CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:),                  &
                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:)))

      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
      i0=nbp_lev-locate(pintou(:),ptropou(i))+1
      ll=check_ozone(v3(i, 1:i0     ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &
                     5.E-9,3.0E-7)
!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
      ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', &
                     5.E-8,1.5E-5)
!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)

    END DO
  END IF

CONTAINS


!-------------------------------------------------------------------------------
!
SUBROUTINE update_fields()
!
!-------------------------------------------------------------------------------
  IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
    CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
    irec=MIN(INT(julien)+1,ntim_in)     !--- irec is just the julian day number
    !--- MIN -> Security in the unlikely case of roundup errors.
    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")
    END IF
  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
    IF(lfirst) irec=locate(time_in,julien) !--- Need to locate surrounding times
    IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
    CALL msg(.TRUE.,'Refreshing adjacent Ozone forcing fields.',sub)
    IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
      WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update (step 1): '&
      //'reading record ',irec,' (time ',time_in(irec),')'
      CALL msg(.TRUE.,str,sub)
      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")
      END IF
    ELSE                                !=== SHIFT FIELDS
      irec=irec+1
      WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update: shifting'&
      //' current next one (',irec,', time ',time_in(irec),')'
      CALL msg(.TRUE.,str,sub)
      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
      END IF
    END IF
    irec=irec+1
    WRITE(str,'(a,i3,a,f12.8,a)')'Next available field update: reading record'&
    ,irec,' (time ',time_in(irec),')'
    CALL msg(.TRUE.,str,sub)
    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")
    END IF
    irec=irec-1
  END IF

END SUBROUTINE update_fields
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE get_2Dfield(v,var)
!
!-------------------------------------------------------------------------------
! Purpose: Shortcut to get the "irec"th record of the surface field named "var"
!          from the input file.
! Remark: In case the field is zonal, it is duplicated along first dimension.
!-------------------------------------------------------------------------------
! Arguments:
  REAL,             INTENT(INOUT) :: v(:,:)
  CHARACTER(LEN=*), INTENT(IN)    :: var
!-------------------------------------------------------------------------------
  CALL NF95_INQ_VARID(fID, TRIM(var), vID)
  CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
  IF(n_dim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
  IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
  CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)

  !--- Flip latitudes: ascending in input file, descending in "rlatu".
  IF(n_dim==3) THEN
    v(1,:) = v(1,nbp_lat:1:-1)
    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
  ELSE
    v(:,:) = v(:,nbp_lat:1:-1)
  END IF

END SUBROUTINE get_2Dfield
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE get_3Dfields(v)
!
!-------------------------------------------------------------------------------
! Purpose: Shortcut to get the "irec"th record of the 3D fields named "nam"
!          from the input file. Fields are stacked on fourth dimension.
! Remark: In case the field is zonal, it is duplicated along first dimension.
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(INOUT) :: v(:,:,:,:)
!-------------------------------------------------------------------------------
  DO i=1,SIZE(nam)
    CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
    IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
    IF(n_dim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
    CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
  END DO

  !--- Flip latitudes: ascending in input file, descending in "rlatu".
  IF(n_dim==3) THEN
    v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
  ELSE
    v(:,:,:,:) = v(:,nbp_lat:1:-1,:,:)
  END IF

END SUBROUTINE get_3Dfields
!
!-------------------------------------------------------------------------------



!-------------------------------------------------------------------------------
!
FUNCTION get_SigTrop(ih,it)
!
!-------------------------------------------------------------------------------
! Arguments:
  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)

END FUNCTION get_SigTrop
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop)
!
!-------------------------------------------------------------------------------
! Purpose: Determine the tropopause using chemical criterium, ie the pressure
!          at which the ozone concentration reaches a certain level.
!-------------------------------------------------------------------------------
! Remarks:
! * Input field is upside down (increasing pressure // increasing vertical idx)
!   The sweep is done from top to ground, starting at the 50hPa layer (idx it0),
!   where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached:
!   the O3 profile in this zone is decreasing with pressure.
! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or
!   determined using an empirical law roughly derived from ... & al.
!-------------------------------------------------------------------------------
! Arguments:
  REAL ::                   chem_tropopause    !--- Pressure at tropopause
  INTEGER,        INTENT(IN)  :: ih            !--- Horizontal index
  INTEGER,        INTENT(OUT) :: it            !--- Index of tropopause layer
  INTEGER,        INTENT(IN)  :: it0           !--- Idx: higher than tropopause
  REAL,           INTENT(IN)  :: pint(:)       !--- Cells-interf Pr, increasing
  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)
  REAL, OPTIONAL, INTENT(IN)  :: pcen(:)       !--- Cells-center Pr, increasing
  REAL, OPTIONAL, INTENT(IN)  :: o3trop        !--- Ozone at tropopause
!-------------------------------------------------------------------------------
! Local variables:
  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
!-------------------------------------------------------------------------------
  !--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN
  IF(PRESENT(o3trop)) THEN                     !=== THRESHOLD FROM ARGUMENTS
    o3t=o3trop
  ELSE                                         !=== EMPIRICAL LAW
    coef = TANH(lat_in(ih)/co3(3))             !--- Latitude  modulation
    coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Seasonal  modulation
    o3t  = 1.E-9 * (co3(1)+co3(2)*coef)        !--- Value from parametrization
  END IF

  !--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)
  it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
  al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
  IF(PRESENT(pcen)) THEN
    chem_tropopause =       pcen(it)**(1.-al) * pcen(it+1)**al
  ELSE
    chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) )
  END IF
  it = locate(pint(:), chem_tropopause)        !--- pint(it)<ptrop<pint(it+1)

END FUNCTION chem_tropopause
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax)
!
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  LOGICAL                      :: check_ozone      !--- .T. => some wrong values
  REAL,             INTENT(IN) :: o3col(:), lon, lat
  INTEGER,          INTENT(IN) :: ilev0
  CHARACTER(LEN=*), INTENT(IN) :: layer
  REAL, OPTIONAL,   INTENT(IN) :: vmin
  REAL, OPTIONAL,   INTENT(IN) :: vmax
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: k
  LOGICAL :: lmin, lmax
  REAL    :: fac
  CHARACTER(LEN=6) :: unt
!-------------------------------------------------------------------------------
  !--- NOTHING TO DO
  lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
  lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
  check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN

  !--- SOME TOO LOW VALUES FOUND
  IF(lmin) THEN
    CALL unitc(vmin,unt,fac)
    DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE
      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', &
        fac*vmin,unt//' in '//TRIM(layer)
    END DO
  END IF

  !--- SOME TOO HIGH VALUES FOUND
  IF(lmax) THEN
    CALL unitc(vmax,unt,fac)
    DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE
      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', &
        fac*vmax,unt//' in '//TRIM(layer)
    END DO
  END IF

END FUNCTION check_ozone
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE unitc(val,unt,fac)
!
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  REAL,             INTENT(IN)  :: val
  CHARACTER(LEN=6), INTENT(OUT) :: unt
  REAL,             INTENT(OUT) :: fac
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: ndgt
!-------------------------------------------------------------------------------
  ndgt=3*FLOOR(LOG10(val)/3.)
  SELECT CASE(ndgt)
    CASE(-6);     unt=' ppmV '; fac=1.E6
    CASE(-9);     unt=' ppbV '; fac=1.E9
    CASE DEFAULT; unt=' vmr  '; fac=1.0
  END SELECT

END SUBROUTINE unitc
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE msg(ll,str,sub)
!
!-------------------------------------------------------------------------------
  USE print_control_mod, ONLY: lunout
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  LOGICAL,                    INTENT(IN) :: ll
  CHARACTER(LEN=*),           INTENT(IN) :: str
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
!-------------------------------------------------------------------------------
  IF(.NOT.ll) RETURN
  IF(PRESENT(sub)) THEN
    WRITE(lunout,*)TRIM(sub)//': '//TRIM(str)
  ELSE
    WRITE(lunout,*)TRIM(str)
  END IF

END SUBROUTINE msg
!
!-------------------------------------------------------------------------------

END SUBROUTINE regr_pr_time_av
!
!-------------------------------------------------------------------------------

END MODULE regr_pr_time_av_m