! $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

  USE print_control_mod,  ONLY: prt_level

  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.

!     - "Ptrp_ou": target grid (LMDZ) tropopause pressure.

!   so that the tropopause is shifted to the position of the one of LMDZ.

!  Note that the ozone quantity conservation is not ensured.

!-------------------------------------------------------------------------------

! 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 centers.

!  * Regridding in pressure can be:

!    - Ploc=='I': pressures provided at Interfaces    => conservative 2nd order

!         OK for ozone coefficients regridding in Cariolle routines.

!    - Ploc=='C': pressures provides at cells Centers => classical linear

!         OK for ozone vertical regridding, especially when torpopause

!      adjustment is activated, to avoid "strairs shape effect" on profiles.

!  * All the fields are regridded as a single multi-dimensional array, so it

!    saves CPU time to call this procedure once for several NetCDF variables

!    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, constant or

!    determined using an empirical three parameters law:

!         o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3)

!       => co1 and co2 are in ppbV, and co3 in degrees.

!       => 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-name "rlonv(i)-rlatu(j)-nam(l)",

!    pressure level/interval (Ploc=="C"/"I") "pcen_in(k)"/"pcen_in(k:k+1)]".

!  * In the 2D file case, the values are the same for all longitudes.

!  * The tropopause correction works like this: the input fields (file) are

!  interpolated on output (LMDZ) pressure field, which is streched using a power

!  law in a limited zone made of 2 layers:

!    1) between lower bound (lower than lowest tropopause) and LMDZ tropopause

!    2) between LMDZ tropopause and upper bound (higher thzn highest tropopause)

!  The stretching function has the following form:

!        Sigma_str = Sigma^(1+alpha*phi(Sigma)), where:

!   * alpha=LOG(SigT_in/SigT_ou)/LOG(SigT_ou)

!     This value shifts the file tropopause to the height of the one of LMDZ.

!   * phi is quasi-linear (sections of 1/log function) in the adjacent intervals:

!       - from 0 to 1 in [Sig_top,SigT_ou]

!       - from 1 to 0 in [SigT_ou,Sig_bot]

!  This quasi-triangular localization function ponderates alpha-law from one near

!  the tropopause to zero each side apart.

!

! * The following fields are on the global "dynamics" grid, as read from files:

  REAL,    SAVE, ALLOCATABLE :: v1 (:,:,:,:)       !--- Current  time ozone fields

! v1: Field read/interpol at time "julien" on the global "dynamics" horiz. grid.

  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous time ozone fields

  REAL,    SAVE, ALLOCATABLE :: v1p(:,:,:,:)       !--- Next     time ozone fields

  REAL,    SAVE, ALLOCATABLE :: pgm(:,:), pgp(:,:) !--- Ground     pressure

  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure

  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio

  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 :: lPrSfile                        !--- Surface pressure flag

  LOGICAL, SAVE :: lPrTfile                        !--- Tropopause pressure flag

  LOGICAL, SAVE :: lO3Tfile                        !--- Tropopause ozone flag

  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag

!$OMP THREADPRIVATE(lfirst)

  REAL,    PARAMETER :: pTropUp=9.E+3 !--- Value  <  tropopause pressure (Pa)

  REAL,    PARAMETER :: gamm   =0.4   !--- Max. stretched layer sigma ratio

  REAL,    PARAMETER :: rho    =1.4   !--- Max tropopauses sigma ratio

  REAL,    PARAMETER :: o3t0   =1.E-7 !--- Nominal O3 vmr at tropopause

  LOGICAL, PARAMETER :: lO3Tpara=.FALSE. !--- Parametrized O3 vmr at tropopause

  LOGICAL, PARAMETER :: ldebug=.FALSE.!--- Force writefield_phy multiple outputs

  REAL, PARAMETER :: ChemPTrMin=9.E+3 !--- Thresholds for minimum and maximum

  REAL, PARAMETER :: ChemPTrMax=4.E+4 !    chemical  tropopause pressure (Pa).

  REAL, PARAMETER :: DynPTrMin =8.E+3 !--- Thresholds for minimum and maximum

  REAL, PARAMETER :: DynPTrMax =4.E+4 !    dynamical tropopause pressure (Pa).

CONTAINS

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

  USE dimphy,         ONLY: klon

  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, &

                            NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION, nf95_get_var

  USE netcdf,         ONLY: NF90_INQ_VARID, 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: nlon=>nbp_lon, nlat=>nbp_lat, nlev_ou=>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

  CHARACTER(LEN=1),  INTENT(IN) :: Ploc       !--- Pressures locations

  !--- File/LMDZ (resp. decreasing & increasing order) pressure, Pa

  !    At cells centers or interfaces depending on "Ploc" keyword (C/I)

  REAL,    INTENT(IN)  :: Pre_in(:)           !--- in:  file      (nlev_in[+1])

  REAL,    INTENT(IN)  :: Pre_ou(:,:)         !--- out: LMDZ (klon,nlev_ou[+1])

  REAL,    INTENT(OUT) :: v3(:,:,:)           !--- Regr. fld (klon,nlev_ou,n_var)

  REAL,    INTENT(IN), OPTIONAL :: Pgnd_ou(:) !--- LMDZ ground pressure   (klon)

  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days

                                              !    since Jan 1 of current year

  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- File longitudes vector (klon)

  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- File latitudes  vector (klon)

  REAL,    INTENT(IN), OPTIONAL :: Ptrp_ou(:) !--- LMDZ 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, iout, nn

  INTEGER :: i, nlev_in, n_dim, itop, itrp, i0

  LOGICAL :: lAdjTro                          !--- Need to adjust tropopause

  REAL    :: y_frac                           !--- Elapsed year fraction

  REAL    :: alpha, beta, al                  !--- For stretching/interpolation

  REAL    :: SigT_in, SigT_ou                 !--- Input and output tropopauses

  REAL    :: Sig_bot, Sig_top                 !--- Bounds of quasi-hat function

  REAL    :: Sig_bo0, Sig_to0                 !--- Lower/upper tropopauses

! * The following fields are scattered to the partial "physics" horizontal grid.

  REAL, POINTER :: v2(:,:,:)                  !--- Current  time ozone fields

!     In the 2D file case, the values are the same for all longitudes.

!     "v2(i, k, l)" is at longitude-latitude "xlon(i)-xlat(i)" and name "nam(l)"

! Both are:          * if Ploc=='I' in pressure interval "press_in_edg(k:k+1)"

!                    * if Ploc=='C' at pressure          "press_in_cen(k)"

  REAL, TARGET :: &

  LOGICAL :: ll

!--- For debug

  REAL, DIMENSION(klon)             :: Ptrop_in, Ptrop_ef

  REAL, DIMENSION(klon)             :: dzStrain, dzStrain0

!-------------------------------------------------------------------------------

  END IF

  !$OMP MASTER

    !=== CHECK WHICH FIELDS ARE AVAILABLE IN THE INPUT FILE

      CALL NF95_INQ_DIMID(fID,"time",vID)

      END IF

      !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)

      ELSE IF(lO3Tpara) THEN; str='    PARAMETRIZED O3 concentration'

    END IF

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

    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES

    END IF

  ELSE

  END IF

  !$OMP END MASTER

  !$OMP BARRIER

  END IF

  ELSE

  ENDIF

  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press

  ELSE

  END IF

  ENDIF

  ENDIF

  !--- No ground pressure in input file => choose it to be the one of LMDZ

  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing

    END DO

!-------------------------------------------------------------------------------

  ELSE                        !--- REGRID IN PRESSURE ; TROPOPAUSE ADJUSTMENT

!-------------------------------------------------------------------------------

    !--- OUTPUT SIGMA LEVELS

      !--- INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS IN CURRENT COLUMN

      !--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered

      ! to keep tropopause pressure realistic ; high values are usually due to

      ! ozone hole fooling the crude chemical tropopause detection algorithm.

      !--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the

      ! dynamical tropopause (especially in filaments) are conterbalanced with

      ! a filter ensuring it stays within a certain distance around input (file)

      ! tropopause, hence avoiding avoid a too thick stretched region ; a final

      ! extra-safety filter keeps the tropopause pressure value realistic.

      !--- POWER LAW COEFFICIENT FOR TROPOPAUSES MATCHING

      !--- DETERMINE STRETCHING DOMAIN UPPER AND LOWER BOUNDS

      END IF

      END IF

      !--- STRETCHING POWER LAW LOCALIZATION FUNCTION:

      !    0 in [0,Sig_top]    0->1 in [Sig_top,SigT_ou]

      !    0 in [Sig_bot,1]    1->0 in [SigT_ou, Sig_bot]

      phi(itop+1:iout) = (1.-LOG(Sig_top)/LOG(Sig_ou(itop+1:iout)))&

      phi(iout+1:ibot) = (1.-LOG(Sig_bot)/LOG(Sig_ou(iout+1:ibot)))&

      !--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)

      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS

      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.

      ll=check_ozone(v3(i, 1:i0-1   ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &

!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)

      ll=check_ozone(v3(i,i0:nlev_ou,1),lon_in(i),lat_in(i),i0,'stratosphere', &

!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)

        dzStrain0(i) = SIGN(7.*LOG(Sig_bo0/Sig_to0),SigT_in-SigT_ou)

        dzStrain (i) = SIGN(7.*LOG(Sig_bot/Sig_top),SigT_in-SigT_ou)

        Ptrop_in (i) = SigT_in*Pgnd_in(i)

        Pstrn_ou(i,:)= Pstr_ou

        phii(i,:)    = phi(:)

        Ptrop_ef(i)  = PTrop_chem(i, itrp, locate(Pres_ou(:),PTropUp),    &

                             Pres_ou(:), v3(:,nlev_ou:1:-1,1),o3trop=o3t0)

      END IF

    END DO

  END IF

  IF(ldebug.AND.lAdjTro) THEN

    CALL writefield_phy('PreSt_ou' ,Pstrn_ou,SIZE(Pre_ou,2)) !--- Strained Pres

    CALL writefield_phy('dzStrain' ,dzStrain ,1)     !--- Strained thickness

    CALL writefield_phy('dzStrain0',dzStrain0,1)     !--- Tropopauses distance

    CALL writefield_phy('phi',phii,nlev_ou)          !--- Localization function

    !--- Tropopauses pressures:

    CALL writefield_phy('PreTr_in',Ptrop_in,1)       !--- Input and effective

    CALL writefield_phy('PreTr_ou',Ptrop_ou,1)       !--- LMDz dyn tropopause

    CALL writefield_phy('PreTr_ef',Ptrop_ef,1)       !--- Effective chem tropop

  END IF

    CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field

    CALL writefield_phy('Ozone_ou',v3(:,:,1),nlev_ou)!--- Output ozone field

    CALL writefield_phy('Pres_ou' ,Pre_ou,SIZE(Pre_ou,2))!--- LMDZ Pressure

  END IF

CONTAINS

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

    !--- MIN -> Security in the unlikely case of roundup errors.

    END IF

  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS

      WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update (step 1): '&

      //'reading record ',irec,' (time ',time_in(irec),')'

      END IF

    ELSE                                !=== SHIFT FIELDS

      WRITE(str,'(a,i3,a,f12.8,a)')'Previous available field update: shifting'&

      //' current next one (',irec,', time ',time_in(irec),')'

      END IF

    END IF

    END IF

  END IF

END SUBROUTINE update_fields

!

!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

! 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_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)

  !--- Flip latitudes: ascending in input file, descending in "rlatu".

  ELSE

  END IF

!

!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

! 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(:,:,:,:)

!-------------------------------------------------------------------------------

    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)

  END DO

  !--- Flip latitudes: ascending in input file, descending in "rlatu".

  ELSE

  END IF

!

!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

! Arguments:

  REAL                 :: out

  INTEGER, INTENT(IN)  :: ih

  INTEGER, INTENT(OUT) :: it

!-------------------------------------------------------------------------------

  !--- Pressure at tropopause read from the forcing file

  !--- Chemical tropopause definition based on a particular threshold

  ELSE IF(lO3Tpara) THEN; out=PTrop_chem(ih,it,itrp0,Pre_in,v2(:,:,1))

!

!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

! 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                        :: out           !--- Pressure at tropopause

  INTEGER,        INTENT(IN)  :: ih            !--- Horizontal index

  INTEGER,        INTENT(OUT) :: it            !--- Index of tropopause layer

  INTEGER,        INTENT(IN)  :: it0           !--- Idx: higher than tropopause

  REAL,           INTENT(IN)  :: pres(:)       !--- Pressure profile, increasing

  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)

  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

  ELSE                                         !=== EMPIRICAL LAW

  END IF

  !--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)

!

!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

  IMPLICIT NONE

!-------------------------------------------------------------------------------

! Arguments:

  LOGICAL                      :: out          !--- .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

  !--- SOME TOO LOW VALUES FOUND

      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//' < ', &

    END DO

  END IF

  !--- SOME TOO HIGH VALUES FOUND

      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//' > ', &

    END DO

  END IF

END FUNCTION check_ozone

!

!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

  IMPLICIT NONE

!-------------------------------------------------------------------------------

! Arguments:

  REAL,             INTENT(IN)  :: val

  CHARACTER(LEN=6), INTENT(OUT) :: unt

  REAL,             INTENT(OUT) :: fac

!-------------------------------------------------------------------------------

! Local variables:

  INTEGER :: ndgt

!-------------------------------------------------------------------------------

  END SELECT

!

!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------

!

!

!-------------------------------------------------------------------------------

  USE print_control_mod, ONLY: lunout

  IMPLICIT NONE

!-------------------------------------------------------------------------------

! Arguments:

  LOGICAL,                    INTENT(IN) :: ll

  CHARACTER(LEN=*),           INTENT(IN) :: str

  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub

!-------------------------------------------------------------------------------

  ELSE

  END IF

END SUBROUTINE msg

!

!-------------------------------------------------------------------------------

END SUBROUTINE regr_pr_time_av

!

!-------------------------------------------------------------------------------

END MODULE regr_pr_time_av_m