Changeset 2968 for LMDZ5


Ignore:
Timestamp:
Jul 25, 2017, 6:03:55 PM (7 years ago)
Author:
dcugnet
Message:

Improved method for ozone forcing: the ozone field from the forcing file is
vertically interpolated on the LMDZ pressure field that is stretched, so
that the resulting tropopause is close from the one of LMDZ.
(need to activate key adjust_tropopause=y in run.def)

The stretching is done in a limited region in each column, defined by a
trapezoidal function:

  • flat between LMDZ and file tropopauses (full stretching)
  • quasi-linear in the 2 surrounding regions (from full to null stretching)

This way, most of the ozone profile is unaffected by the stretching.

The total ozone column is not stricly conserved.

The transition regions are thick enough to ensure that the stretched
pressure profile is growing with decreasing height. Some additional margin
helps to avoid a too steep variation in the stretched region.

In case the distance separating the chemical tropopause of the forcing
fields (built from monthly means) from the LMDZ dynamical tropopause
(extremely variable, both in space and time) is too large, it is under-
evaluated to avoid injecting ozone at low troposphere levels (filaments
of stratospheric air entering the troposphere are not accurately enough
detected to allow that to be done reliably).

In the dynamical tropopause detection routine, the PV and the potential
temperature are moderately smoothed to avoid the algorithm to be fooled
by local variations.

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

Legend:

Unmodified
Added
Removed

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

    r2952 r2968  
    195195       beta_prec,  &
    196196       rneb,  &
    197        zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
    198        pr_tropopause
     197       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
    199198       !
    200199    USE phys_state_var_mod ! Variables sauvegardees de la physique
    … …  
    20152014          IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
    20162015          IF(adjust_tropopause) THEN
    2017              CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot, &
    2018                   pr_tropopause)
    20192016             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
    20202017                      ro3i,  press_edg_climoz,  paprs,   wo, time_climoz,    &
    2021                            latitude_deg,  press_cen_climoz,  pr_tropopause)
     2018                      longitude_deg,   latitude_deg,    press_cen_climoz,    &
     2019                      dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
    20222020          ELSE
    20232021             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &

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

    r2820 r2968  
    22MODULE regr_pr_time_av_m
    33
    4   USE print_control_mod, ONLY: lunout
    54  USE write_field_phy
    65  USE mod_phys_lmdz_transfert_para, ONLY: bcast
     6  USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank
    77  IMPLICIT NONE
    88
    … …  
    1414!  * in all the threads of all the processes, regrid the fields in pressure
    1515!    to the LMDZ vertical grid.
    16 !  * input files fields are stretched in the viscinity (+/- 5 kms) of the mean
    17 !    tropopause (geometrical mean of LMDZ and input fields tropopauses) to force
    18 !    the resulting field tropopause height to the one of LMDZ. To switch this
    19 !    feature on, the following arguments must be present:
     16!  * the forcing fields are stretched if the following arguments are present:
    2017!     - "lat_in":   input file latitudes.
    2118!     - "pcen_in":  input file cells center pressure levels.
    2219!     - "ptrop_ou": target grid (LMDZ) tropopause pressure.
    23 !    Note that the ozone quantity conservation is not ensured.
     20!   so that the tropopause is shifted to the position of the one of LMDZ.
     21!  Note that the ozone quantity conservation is not ensured.
    2422!-------------------------------------------------------------------------------
    2523! Initial routine: regr_pr_av_m module (L. Guez).
    … …  
    2725!    - time interpolation
    2826!    - 3D compliant
    29 !    - vertical stretching of the field to allow tropopause and ground pressure
    30 !      matching between input field and current lmdz field.
     27!    - vertical stretching of the field to allow tropopause matching between
     28!    input field and current lmdz field.
    3129!-------------------------------------------------------------------------------
    3230! Remarks:
    … …  
    4341!  * The input file pressure at tropopause can be (in decreasing priority):
    4442!    1) read from the file if "tropopause_air_pressure" field is available.
    45 !    2) computed from the input file ozone field using:
    46 !       - o3 concentration at tropopause if "tro3_at_tropopause" is available.
    47 !       - a default value (100ppbv) if not.
     43!    2) computed using "tro3" and "tro3_at_tropopause' (if available).
     44!    3) computed using "tro3" and a fixed threshold otherwise, determined using
     45!    an empirical three parameters law:
     46!         o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/co3)
     47!       => co1 and co2 are in ppbV, and co3 in degrees.
     48!       => co3 allow a smooth transition between north and south hemispheres.
    4849!  * If available, the field "ps" (input file pressure at surface) is used.
    4950!    If not, the current LMDZ ground pressure is taken instead.
    50 !  * The O3 threshold for tropopause is defined using 3 coefficients:
    51 !    o3t(ppbV)=co1+co2*SIN(PI*(month-2)/6)*TANH(lat_deg/c)
    52 !     => co1 and co2 are in ppmV, and co3 in degrees.
    53 !     => co3 allow a smooth transition between north and south hemispheres.
    54 !-------------------------------------------------------------------------------
    55 ! * Fields with suffix "m"/"p" are at the closest records earlier/later than the
    56 !   middle of the julian day "julien", on the global "dynamics" horizontal grid.
    57 ! * Fields(i,j,k,l) are at longitude-latitude "rlonv(i)-rlatu(j)", pressure
    58 !   interval "pint_in(k:k+1)]" and variable "nam(l)".
    59 ! * In the 2D file case, the values are the same for all longitudes.
     51!  * Fields with suffix "m"/"p" are at the closest records earlier/later than
     52!  the mid-julian day "julien", on the global "dynamics" horizontal grid.
     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)".
     55!  * In the 2D file case, the values are the same for all longitudes.
     56!  * The tropopause correction works like this: the input fields (file) are
     57!  interpolated on output (LMDZ) pressure field, which is streched using a power
     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
     61!  The stretching function has the following form:
     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.
     70!
    6071! * The following fields are on the global "dynamics" grid, as read from files:
    6172  REAL,    SAVE, ALLOCATABLE :: v1m(:,:,:,:)       !--- Previous ozone fields
    … …  
    6475  REAL,    SAVE, ALLOCATABLE :: ptm(:,:), ptp(:,:) !--- Tropopause pressure
    6576  REAL,    SAVE, ALLOCATABLE :: otm(:,:), otp(:,:) !--- Tropopause o3 mix. ratio
     77  INTEGER, SAVE :: ntim_in                         !--- Records nb in input file
     78  INTEGER, SAVE :: itrp0                           !--- idx above chem tropop.
    6679  INTEGER, SAVE :: irec                            !--- Current time index
    6780!      * for daily   input files: current julian day number
    … …  
    7386  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
    7487!$OMP THREADPRIVATE(lfirst)
    75   REAL, PARAMETER :: co3(3)=[91.,28.,20.]          !--- Coeffs for o3 threshold
    76   REAL, PARAMETER :: prtrop=5.E+3  !--- Value lower than the tropopause pressure
    77   REAL, PARAMETER :: delta =5.     !--- Dist. around tropopause for strain (kms)
     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 paramztrized O3 vmr at tropopause
    7893
    7994CONTAINS
    … …  
    8196!-------------------------------------------------------------------------------
    8297!
    83 SUBROUTINE regr_pr_time_av(fID, nam, julien, pint_in, pint_ou, v3, time_in,    &
    84                                      lat_in, pcen_in, ptrop_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)
    85100!
    86101!-------------------------------------------------------------------------------
    87102  USE dimphy,         ONLY: klon
    88   USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err
     103  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, handle_err, &
     104                            NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION
    89105  USE netcdf,         ONLY: NF90_INQ_VARID, NF90_GET_VAR, NF90_NOERR
    90106  USE assert_m,       ONLY: assert
    … …  
    105121  REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
    106122  REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces file Pres, Pa, ascending
    107   REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ Pres, Pa (klon,llm+1)
     123  REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ pressure, Pa (g,nbp_lev+1)
    108124  REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
    109125  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
    110126                                              !    since Jan 1 of current year
     127  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- Input/output longitudes vector
    111128  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
    112129  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure
    113   REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pressure (klon)
     130  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pres (klon)
    114131!-------------------------------------------------------------------------------
    115132! Local variables:
    … …  
    117134  include "YOMCST.h"
    118135  CHARACTER(LEN=80)  :: sub
    119   CHARACTER(LEN=320) :: msg
    120   INTEGER :: vID, ncerr, n_var, nlev_in,ntim_in, ndim, i, ibot, itop, itrp,itrp0
     136  CHARACTER(LEN=320) :: str
     137  INTEGER :: vID, ncerr, n_var, ibot, ibo0, nn, itrp
     138  INTEGER :: i, nlev_in, n_dim, itop, ito0, i0
    121139  LOGICAL :: lAdjTro                          !--- Need to adjust tropopause
    122140  REAL    :: y_frac                           !--- Elapsed year fraction
    123   REAL    :: alpha, beta, al                  !--- For strectching/interpolation
    124   REAL    :: SigT_in, SigT_ou, SigT_mn        !--- Tropopause: in/out/mean
    125   REAL    :: SigA_bot, SigA_top               !--- Strained domain bounds
    126   REAL    :: Sig_in (SIZE(pint_in))           !--- Input field sigma levels
    127   REAL    :: phi    (SIZE(pint_in))           !--- Stretching exponent anomaly
    128   REAl    :: Pint_st(SIZE(pint_in))           !--- Stretched pressure levels
     141  REAL    :: alpha, beta, al                  !--- For stretching/interpolation
     142  REAL    :: SigT_in, SigT_ou                 !--- Input and output tropopauses
     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
    129150  REAL    :: v1(nbp_lon,nbp_lat,SIZE(pint_in)-1,SIZE(nam))
    130151  REAL    :: v2(klon,           SIZE(pint_in)-1,SIZE(nam))
    … …  
    134155!     "v2(i,    k, l)" is at longitude-latitude  "xlon(i)-xlat(i)".
    135156! Both are for pressure interval "press_in_edg(k:k+1)]" and variable "nam(l)"
    136   REAL, DIMENSION(nbp_lon, nbp_lat) ::  ps1, pt1, ot1
    137   REAL, DIMENSION(klon)             ::  ps2, pt2, ot2
     157  REAL, DIMENSION(nbp_lon, nbp_lat)   :: ps1, pt1, ot1
     158  REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
     159  LOGICAL :: ll
    138160!-------------------------------------------------------------------------------
    139161  sub="regr_pr_time_av"
    140   nlev_in=SIZE(pint_in)-1; ntim_in=SIZE(time_in)
    141   CALL assert(SIZE(v3,1)==klon,                 TRIM(sub)//" v3 klon")
    142   CALL assert(SIZE(v3,2)==nbp_lev,              TRIM(sub)//" v3 nbp_lev")
    143   n_var = assert_eq(SIZE(nam), SIZE(v3,3),      TRIM(sub)//" v3 n_var")
    144   CALL assert(SHAPE(pint_ou)==[klon, nbp_lev+1],TRIM(sub)//" pint_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")
     165  n_var = assert_eq(SIZE(nam),SIZE(v3,3),TRIM(sub)//" v3 n_var")
     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")
    145169  IF(PRESENT(lat_in))   CALL assert(SIZE(lat_in  )==klon,TRIM(sub)//" lat_in klon")
    146170  IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon")
    147   IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in)==nlev_in,TRIM(sub)//" pcen_in")
     171  IF(PRESENT(pcen_in))  CALL assert(SIZE(pcen_in )==nlev_in,TRIM(sub)//" pcen_in")
    148172  lAdjTro=PRESENT(ptrop_ou)
    149173  IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) &
    … …  
    158182      lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR
    159183      lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==NF90_NOERR
    160       linterp=PRESENT(time_in); IF(linterp) linterp=ntim_in==14
     184      CALL NF95_INQ_DIMID(fID,"time",vID)
     185      CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
     186      linterp=PRESENT(time_in).AND.ntim_in==14
    161187      IF(linterp) THEN
    162188        ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var))
    … …  
    167193      END IF
    168194      !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa)
    169       IF(lAdjTro) itrp0=locate(pcen_in,prtrop)
    170       IF(lPrSurf) WRITE(lunout,*)TRIM(sub)//': Using GROUND PRESSURE from input O3 forcing file.'
    171       IF(linterp) WRITE(lunout,*)TRIM(sub)//': Monthly O3 files => ONLINE TIME INTERPOLATION.'
    172       IF(lAdjTro) THEN
    173         WRITE(lunout,*)TRIM(sub)//': o3 forcing file tropopause location uses:'
    174         IF(lPrTrop) THEN
    175           WRITE(lunout,*)'    PRESSURE AT TROPOPAUSE from file.'
    176         ELSE IF(lO3Trop) THEN
    177           WRITE(lunout,*)'    O3 CONCENTRATION AT TROPOPAUSE from file.'
    178         ELSE
    179           WRITE(lunout,*)'    PARAMETRIZED O3 concentration at tropopause.'
    180         END IF
     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
     204        CALL msg(lAdjTro,'    PARAMETRIZED O3 concentration at tropopause.')
    181205      END IF
    182206    END IF
    … …  
    185209    CALL update_fields()
    186210
    187 
    188211    !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES
    189212    IF(linterp) THEN
    190       WRITE(lunout,'(3(a,f7.3))')TRIM(sub)//': Interpolating O3 at julian day '&
    191       &,julien,' from fields at times ',time_in(irec),' and ', time_in(irec+1)
     213      WRITE(str,'(3(a,f7.3))')'Interpolating O3 at julian day ',julien,' from '&
     214      &//'fields at times ',time_in(irec),' and ', time_in(irec+1)
     215      CALL msg(.TRUE.,str,sub)
    192216      al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec))
    193217      v1=al*v1m+(1.-al)*v1p
    … …  
    199223  !$OMP END MASTER
    200224  IF(lfirst) THEN
    201     lfirst=.FALSE. ;     CALL bcast(lfirst) ; CALL bcast(lPrSurf)
    202     CALL bcast(lPrTrop); CALL bcast(lO3Trop); CALL bcast(linterp)
    203     IF(lAdjTro) CALL bcast(itrp0)
     225    lfirst=.FALSE.;      CALL bcast(lfirst)
     226    IF(lAdjTro)          CALL bcast(itrp0)
     227    CALL bcast(lPrTrop); CALL bcast(lPrSurf)
     228    CALL bcast(lO3Trop); CALL bcast(linterp)
    204229  END IF
    205230  CALL scatter2d(v1,v2)
    206231  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
    207   IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=Pint_ou(:,1); END IF
     232  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF
    208233  IF(lPrTrop) CALL scatter2d(pt1,pt2)
    209234  IF(lO3Trop) CALL scatter2d(ot1,ot2)
    … …  
    212237  IF(.NOT.lAdjTro) THEN
    213238    DO i=1,klon
    214       CALL regr_conserv(1,v2(i,:,:) , Pint_in , Pint_ou(i,nbp_lev+1:1:-1),    &
    215                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in))
     239      pintou = pint_ou(i,nbp_lev+1:1:-1)
     240      CALL regr_conserv(1,v2(i,:,:), pint_in(:), pintou(:),                    &
     241                          v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:), pint_in(:)))
    216242    END DO
    217243  ELSE
    218244    y_frac=(REAL(days_elapsed)+jH_cur)/year_len
     245
     246    !--- OUTPUT SIGMA LEVELS
    219247    DO i=1,klon
    220       SigT_in = get_SigTrop(i,itrp)        !--- input  (file) tropopause
    221       SigT_ou = ptrop_ou(i)/pint_ou(i,1)   !--- output (lmdz) tropopause
    222       SigT_mn = SQRT(SigT_in*SigT_ou)      !--- mean tropopause>strained domain
    223 
    224       !--- DEFINE THE VERTICAL LAYER IN WHICH THE PROFILE IS STRECHED
    225       beta=EXP(delta/scaleheight); Sig_in(:) = [pint_in(1:nlev_in)/ps2(i),1.]
    226       SigA_bot = SigT_mn*beta ; ibot=locate(Sig_in(:),SigA_bot)
    227       SigA_top = SigT_mn/beta ; itop=locate(Sig_in(:),SigA_top)+1
    228 
    229       !--- HAT FUNCTION phi (/=0 in [SigA_bot,SigA_top] only)
     248
     249      !--- LOCAL INPUT/OUTPUT (FILE/LMDZ) SIGMA LEVELS AT INTERFACES
     250      pintou=pint_ou(i,nbp_lev+1:1:-1)            !--- increasing values
     251      Sig_in(:) = [pint_in(1:nlev_in+1)/ps2(i)]   !--- increasing values
     252      Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values
     253
     254      !--- INPUT (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE
     255      SigT_in = get_SigTrop(i,itrp)
     256      SigT_ou = ptrop_ou(i)/ps2(i)
     257
     258      !--- AVOID THE FILAMENTS WHICH WOULD NEED A VERY THICK STRETCHED REGION
     259      IF(SigT_ou>SigT_in*rho) SigT_ou = SigT_in*rho
     260      IF(SigT_ou<SigT_in/rho) SigT_ou = SigT_in/rho
     261      ptropou(i)=SigT_ou*ps2(i)
     262
     263      !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW
     264      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
     265
     266      !--- FULLY STRETCHED LAYER BOUNDS (FILE AND MODEL TROPOPAUSES)
     267      Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot)
     268      Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top)
     269
     270      !--- PARTIALLY STRETCHED LAYER BOUNDS, ENSURING >0 DERIVATIVE
     271      beta = LOG(Sig_top)/LOG(Sig_bot)
     272      Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
     273      Sig_to0 = Sig_top ; IF(alpha>0.) Sig_to0 = Sig_top **  beta
     274
     275      !--- SOME ADDITIONAL MARGIN, PROPORTIONAL TO STRETCHED REGION THICKNESS
     276      !--- gamma<log(Sig_bo0/|alpha|) to keep Sig_bo0<1
     277      Sig_bo0 = MIN(Sig_bo0*EXP( gamm*ABS(alpha)), 0.95+(1.-0.95)*Sig_bo0)
     278      Sig_to0 =     Sig_to0*EXP(-gamm*ABS(alpha))
     279      ibo0 = locate(Sig_ou(:),Sig_bo0)
     280      ito0 = locate(Sig_ou(:),Sig_to0)
     281
     282      !--- FUNCTION FOR STRETCHING LOCALISATION
     283      !    0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1
    230284      phi(:)=0.
    231       phi(itop:itrp)=(Sig_in(itop:itrp)-SigA_top)/(SigT_in-SigA_top)
    232       phi(itrp:ibot)=(SigA_bot-Sig_in(itrp:ibot))/(SigA_bot-SigT_in)
    233 
    234       !--- STRAINED INPUT logP PROFILE ; alpha: MAX. STRETCH. EXPONENT INCREMENT
    235       alpha = LOG(SigT_ou/SigT_in)/LOG(SigT_in)
    236       Pint_st(:) = pint_ou(i,1) * Sig_in(:)**(1+alpha*phi(:))
    237 
    238       !--- REGRID INPUT PROFILE ON STRAINED VERTICAL LEVELS
    239       CALL regr_conserv(1,v2(i,:,:) , Pint_st,  Pint_ou(i,nbp_lev+1:1:-1),     &
    240                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),Pint_st))
     285      phi(itop+1:ibot) =  1.
     286      phi(ito0+1:itop) = (1.-LOG(Sig_to0)/LOG(Sig_ou(ito0+1:itop)))&
     287                            *LOG(Sig_top)/LOG(Sig_top/Sig_to0)
     288      phi(ibot+1:ibo0) = (1.-LOG(Sig_bo0)/LOG(Sig_ou(ibot+1:ibo0)))&
     289                            *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)
     290
     291      !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
     292      pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
     293
     294      !--- REGRID INPUT PROFILE ON STRAINED VERTICAL OUTPUT LEVELS
     295      CALL regr_conserv(1, v2(i,:,:), pint_in(:), pstr_ou(:),                  &
     296                           v3(i,nbp_lev:1:-1,:), slopes(1,v2(i,:,:),pint_in(:)))
     297
     298      !--- CHECK CONCENTRATIONS. strato: 50ppbV-15ppmV ; tropo: 5ppbV-300ppbV.
     299      i0=nbp_lev-locate(pintou(:),ptropou(i))+1
     300      ll=check_ozone(v3(i, 1:i0     ,1),lon_in(i),lat_in(i),1 ,'troposphere',  &
     301                     5.E-9,3.0E-7)
     302!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in troposphere', 1)
     303      ll=check_ozone(v3(i,i0:nbp_lev,1),lon_in(i),lat_in(i),i0,'stratosphere', &
     304                     5.E-8,1.5E-5)
     305!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
     306
    241307    END DO
    242308  END IF
    243309
    244 
    245310CONTAINS
    246311
    … …  
    252317!-------------------------------------------------------------------------------
    253318  IF(.NOT.linterp) THEN                 !=== DAILY FILES: NO TIME INTERPOLATION
    254     WRITE(lunout,*)TRIM(sub)//': Updating Ozone forcing field: read from file.'
     319    CALL msg(.TRUE.,sub,'Updating Ozone forcing field: read from file.')
    255320    irec=MIN(INT(julien)+1,ntim_in)     !--- irec is just the julian day number
    256321    !--- MIN -> Security in the unlikely case of roundup errors.
    … …  
    263328  ELSE                                  !=== MONTHLY FILES: GET 2 NEAREST RECS
    264329    IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN
    265     WRITE(lunout,*)TRIM(sub)//': Refreshing adjacent Ozone forcing fields.'
     330    CALL msg(.TRUE.,sub,'Refreshing adjacent Ozone forcing fields.')
    266331    IF(lfirst) THEN                     !=== READ EARLIEST TIME FIELDS
    267332      irec=locate(time_in,julien)       !--- Need to locate surrounding times
    … …  
    309374!-------------------------------------------------------------------------------
    310375  CALL NF95_INQ_VARID(fID, TRIM(var), vID)
    311   CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=ndim)
    312   IF(ndim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
    313   IF(ndim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
     376  CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
     377  IF(n_dim==2) ncerr=NF90_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
     378  IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(:,:), start=[1,1,irec])
    314379  CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(var),ncerr,fID)
    315380
    316381  !--- Flip latitudes: ascending in input file, descending in "rlatu".
    317   IF(ndim==3) THEN
     382  IF(n_dim==3) THEN
    318383    v(1,:) = v(1,nbp_lat:1:-1)
    319384    v(2:,:)= SPREAD(v(1,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
    … …  
    341406  DO i=1,SIZE(nam)
    342407    CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
    343     CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=ndim)
    344     IF(ndim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
    345     IF(ndim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
     408    CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
     409    IF(n_dim==3) ncerr=NF90_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
     410    IF(n_dim==4) ncerr=NF90_GET_VAR(fID,vID,v(:,:,:,i), start=[1,1,1,irec])
    346411    CALL handle_err(TRIM(sub)//" NF90_GET_VAR "//TRIM(nam(i)),ncerr,fID)
    347412  END DO
    348413
    349414  !--- Flip latitudes: ascending in input file, descending in "rlatu".
    350   IF(ndim==3) THEN
     415  IF(n_dim==3) THEN
    351416    v(1,:,:,:) = v(1,nbp_lat:1:-1,:,:)
    352417    v(2:,:,:,:)= SPREAD(v(1,:,:,:),DIM=1,ncopies=nbp_lon-1)  !--- Duplication
    … …  
    360425
    361426
     427
    362428!-------------------------------------------------------------------------------
    363429!
    … …  
    370436  REAL                 :: get_Sigtrop
    371437!-------------------------------------------------------------------------------
     438  !--- Pressure at tropopause is read in the forcing file
     439  IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
     440    get_SigTrop=pt2(ih)/ps2(ih); RETURN
     441  END IF
     442  !--- Chemical tropopause definition is used using a particular threshold
     443  IF(lO3Trop) THEN                             !--- o3trop KNOWN FROM FILE
     444    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,ot2(ih))
     445  ELSE IF(lo3tp) THEN                          !--- o3trop PARAMETRIZATION
     446    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in)
     447  ELSE                                         !--- o3trop CONSTANT
     448    get_SigTrop=chem_tropopause(ih,it,itrp0,pint_in,v2(:,:,1),pcen_in,o3t0)
     449  END IF
     450  get_SigTrop=get_SigTrop/ps2(ih)
     451
     452END FUNCTION get_SigTrop
     453!
     454!-------------------------------------------------------------------------------
     455
     456
     457!-------------------------------------------------------------------------------
     458!
     459FUNCTION chem_tropopause(ih,it,it0,pint,o3,pcen,o3trop)
     460!
     461!-------------------------------------------------------------------------------
     462! Purpose: Determine the tropopause using chemical criterium, ie the pressure
     463!          at which the ozone concentration reaches a certain level.
     464!-------------------------------------------------------------------------------
     465! Remarks:
     466! * Input field is upside down (increasing pressure // increasing vertical idx)
     467!   The sweep is done from top to ground, starting at the 50hPa layer (idx it0),
     468!   where O3 is about 1,5 ppmV, until the first layer where o3<o3t is reached:
     469!   the O3 profile in this zone is decreasing with pressure.
     470! * Threshold o3t can be given as an input argument ("o3trop", in ppmV) or
     471!   determined using an empirical law roughly derived from ... & al.
     472!-------------------------------------------------------------------------------
     473! Arguments:
     474  REAL ::                   chem_tropopause    !--- Pressure at tropopause
     475  INTEGER,        INTENT(IN)  :: ih            !--- Horizontal index
     476  INTEGER,        INTENT(OUT) :: it            !--- Index of tropopause layer
     477  INTEGER,        INTENT(IN)  :: it0           !--- Idx: higher than tropopause
     478  REAL,           INTENT(IN)  :: pint(:)       !--- Cells-interf Pr, increasing
     479  REAL,           INTENT(IN)  :: o3(:,:)       !--- Ozone field (pptV)
     480  REAL, OPTIONAL, INTENT(IN)  :: pcen(:)       !--- Cells-center Pr, increasing
     481  REAL, OPTIONAL, INTENT(IN)  :: o3trop        !--- Ozone at tropopause
     482!-------------------------------------------------------------------------------
    372483! Local variables:
    373   INTEGER :: ii                                !--- Idx of layer containing o3t
    374484  REAL    :: o3t                               !--- Ozone concent. at tropopause
    375   REAL    :: prt                               !--- Air pressure   at tropopause
    376485  REAL    :: al                                !--- Interpolation coefficient
    377   REAL    :: coef                              !--- Coef: North/South transition
    378 !-------------------------------------------------------------------------------
    379   !--- DETERMINE PRESSURE AT TROPOPAUSE AND DIVIDE IT BY GROUND PRESSURE
    380   IF(lPrTrop) THEN                             !--- PrTrop KNOWN FROM FILE
    381     get_SigTrop=pt2(ih)/ps2(ih)
    382     it=locate(pint_in(:),pt2(ih))
    383   ELSE                                         !--- O3 THRESHOLD
    384     coef = TANH(lat_in(i)/co3(3))              !--- Latitude dependant ponderat.
    385     coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Time-dependant ponderation
     486  REAL    :: coef                              !--- Coeff of latitude modulation
     487  REAL, PARAMETER :: co3(3)=[91.,28.,20.]      !--- Coeff for o3 at tropopause
     488!-------------------------------------------------------------------------------
     489  !--- DETERMINE THE OZONE CONCENTRATION AT TROPOPAUSE IN THE CURRENT COLUMN
     490  IF(PRESENT(o3trop)) THEN                     !=== THRESHOLD FROM ARGUMENTS
     491    o3t=o3trop
     492  ELSE                                         !=== EMPIRICAL LAW
     493    coef = TANH(lat_in(ih)/co3(3))             !--- Latitude  modulation
     494    coef = SIN (2.*RPI*(y_frac-1./6.)) * coef  !--- Seasonal  modulation
    386495    o3t  = 1.E-9 * (co3(1)+co3(2)*coef)        !--- Value from parametrization
    387     IF(lO3Trop) o3t=ot2(ih)                    !--- Value from file
    388     !--- Starts from 50hPa and go down.
    389     it=itrp0; DO WHILE(v2(ih,it+1,1)>=o3t); it=it+1; END DO
    390     al=(v2(ih,it,1)-o3t)/(v2(ih,it,1)-v2(ih,it+1,1))
    391     get_SigTrop = ( pcen_in(it)**(1.-al) * pcen_in(it+1)**al )/ps2(ih)
    392   END IF
    393 
    394 END FUNCTION get_SigTrop
    395 !
    396 !-------------------------------------------------------------------------------
    397 
     496  END IF
     497
     498  !--- FROM 50hPa, GO DOWN UNTIL "it" IS SUCH THAT o3(ih,it)>o3t>o3(ih,it+1)
     499  it=it0; DO WHILE(o3(ih,it+1)>=o3t); it=it+1; END DO
     500  al=(o3(ih,it)-o3t)/(o3(ih,it)-o3(ih,it+1))
     501  IF(PRESENT(pcen)) THEN
     502    chem_tropopause =       pcen(it)**(1.-al) * pcen(it+1)**al
     503  ELSE
     504    chem_tropopause = SQRT( pint(it)**(1.-al) * pint(it+2)**al * pint(it+1) )
     505  END IF
     506  it = locate(pint(:), chem_tropopause)        !--- pint(it)<ptrop<pint(it+1)
     507
     508END FUNCTION chem_tropopause
     509!
     510!-------------------------------------------------------------------------------
     511
     512
     513!-------------------------------------------------------------------------------
     514!
     515FUNCTION check_ozone(o3col, lon, lat, ilev0, layer, vmin, vmax)
     516!
     517!-------------------------------------------------------------------------------
     518  IMPLICIT NONE
     519!-------------------------------------------------------------------------------
     520! Arguments:
     521  LOGICAL                      :: check_ozone      !--- .T. => some wrong values
     522  REAL,             INTENT(IN) :: o3col(:), lon, lat
     523  INTEGER,          INTENT(IN) :: ilev0
     524  CHARACTER(LEN=*), INTENT(IN) :: layer
     525  REAL, OPTIONAL,   INTENT(IN) :: vmin
     526  REAL, OPTIONAL,   INTENT(IN) :: vmax
     527!-------------------------------------------------------------------------------
     528! Local variables:
     529  INTEGER :: k
     530  LOGICAL :: lmin, lmax
     531  REAL    :: fac
     532  CHARACTER(LEN=6) :: unt
     533!-------------------------------------------------------------------------------
     534  !--- NOTHING TO DO
     535  lmin=.FALSE.; IF(PRESENT(vmin)) lmin=COUNT(o3col<vmin)/=0
     536  lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
     537  check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN
     538
     539  !--- SOME TOO LOW VALUES FOUND
     540  IF(lmin) THEN
     541    CALL unitc(vmin,unt,fac)
     542    DO k=1,SIZE(o3col); IF(o3col(k)>vmin) CYCLE
     543      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
     544        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' < ', &
     545        fac*vmin,unt//' in '//TRIM(layer)
     546    END DO
     547  END IF
     548
     549  !--- SOME TOO HIGH VALUES FOUND
     550  IF(lmax) THEN
     551    CALL unitc(vmax,unt,fac)
     552    DO k=1,SIZE(o3col); IF(o3col(k)<vmax) CYCLE
     553      WRITE(*,'(a,2f7.1,i3,a,2(f8.3,a))')'WARNING: inconsistent value in '     &
     554        //TRIM(layer)//': O3(',lon,lat,k+ilev0-1,')=',fac*o3col(k),unt//' > ', &
     555        fac*vmax,unt//' in '//TRIM(layer)
     556    END DO
     557  END IF
     558
     559END FUNCTION check_ozone
     560!
     561!-------------------------------------------------------------------------------
     562
     563
     564!-------------------------------------------------------------------------------
     565!
     566SUBROUTINE unitc(val,unt,fac)
     567!
     568!-------------------------------------------------------------------------------
     569  IMPLICIT NONE
     570!-------------------------------------------------------------------------------
     571! Arguments:
     572  REAL,             INTENT(IN)  :: val
     573  CHARACTER(LEN=6), INTENT(OUT) :: unt
     574  REAL,             INTENT(OUT) :: fac
     575!-------------------------------------------------------------------------------
     576! Local variables:
     577  INTEGER :: ndgt
     578!-------------------------------------------------------------------------------
     579  ndgt=3*FLOOR(LOG10(val)/3.)
     580  SELECT CASE(ndgt)
     581    CASE(-6);     unt=' ppmV '; fac=1.E6
     582    CASE(-9);     unt=' ppbV '; fac=1.E9
     583    CASE DEFAULT; unt=' vmr  '; fac=1.0
     584  END SELECT
     585
     586END SUBROUTINE unitc
     587!
     588!-------------------------------------------------------------------------------
     589
     590
     591!-------------------------------------------------------------------------------
     592!
     593SUBROUTINE msg(ll,str,sub)
     594!
     595!-------------------------------------------------------------------------------
     596  USE print_control_mod, ONLY: lunout
     597  IMPLICIT NONE
     598!-------------------------------------------------------------------------------
     599! Arguments:
     600  LOGICAL,                    INTENT(IN) :: ll
     601  CHARACTER(LEN=*),           INTENT(IN) :: str
     602  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
     603!-------------------------------------------------------------------------------
     604  IF(.NOT.ll) RETURN
     605  IF(PRESENT(sub)) THEN
     606    WRITE(lunout,*)TRIM(sub)//': '//TRIM(str)
     607  ELSE
     608    WRITE(lunout,*)TRIM(str)
     609  END IF
     610
     611END SUBROUTINE msg
     612!
     613!-------------------------------------------------------------------------------
    398614
    399615END SUBROUTINE regr_pr_time_av

  • LMDZ5/trunk/libf/phylmd/tropopause_m.F90

    r2788 r2968  
    22
    33!  USE phys_local_var_mod, ONLY: ptrop => pr_tropopause
     4  IMPLICIT NONE
    45  PRIVATE
    56  PUBLIC :: dyn_tropopause
    … …  
    910!-------------------------------------------------------------------------------
    1011!
    11 SUBROUTINE dyn_tropopause(t,ts,paprs,pplay,presnivs,rot,ptrop,tpot0,pvor0,pmin0,pmax0)
     12FUNCTION dyn_tropopause(t, ts, paprs, pplay, rot, thet0, pvor0)
    1213!
    1314!-------------------------------------------------------------------------------
    14   USE assert_m,      ONLY: assert
    15   USE assert_eq_m,   ONLY: assert_eq
    16   USE comvert_mod,   ONLY: preff
    17   USE dimphy,        ONLY: klon, klev
    18   USE geometry_mod,  ONLY: latitude_deg, longitude_deg
    19   USE interpolation, ONLY: locate
    20   IMPLICIT NONE
     15  USE assert_m,     ONLY: assert
     16  USE assert_eq_m,  ONLY: assert_eq
     17  USE dimphy,       ONLY: klon, klev
     18  USE geometry_mod, ONLY: latitude_deg, longitude_deg
     19  USE write_field_phy
     20  USE vertical_layers_mod, ONLY: aps, bps, preff
     21
    2122!-------------------------------------------------------------------------------
    2223! Arguments:
     24  REAL ::     dyn_tropopause(klon) !--- Pressure at tropopause
    2325  REAL, INTENT(IN)  ::      t(:,:) !--- Cells-centers temperature
    2426  REAL, INTENT(IN)  ::     ts(:)   !--- Surface       temperature
    2527  REAL, INTENT(IN)  ::  paprs(:,:) !--- Cells-edges   pressure
    2628  REAL, INTENT(IN)  ::  pplay(:,:) !--- Cells-centers pressure
    27   REAL, INTENT(IN)  :: presnivs(:) !--- Cells-centers nominal pressure
    2829  REAL, INTENT(IN)  ::    rot(:,:) !--- Cells-centers relative vorticity
    29   REAL, INTENT(OUT) :: ptrop(klon) !--- Tropopause pressure
    30   REAL, INTENT(IN), OPTIONAL :: tpot0, pvor0, pmin0, pmax0
     30  REAL, INTENT(IN), OPTIONAL :: thet0, pvor0
    3131!-------------------------------------------------------------------------------
    3232! Local variables:
    3333  include "YOMCST.h"
    3434  CHARACTER(LEN=80)  :: sub
    35   INTEGER :: it, i, nsrf, k, nh, kmin, kmax
    36   REAL    :: p1, p2, tpo0, pvo0, pmn0, pmx0, al
    37   REAL, DIMENSION(klon,klev)   :: t_edg, t_pot
    38   REAL, DIMENSION(klon,klev-1) :: pvort
     35  INTEGER :: i, k, kb, kt, kp, ib, ie, nw
     36  REAL    :: al, th0, pv0
     37  REAL,    DIMENSION(klon,klev) :: tpot_cen, tpot_edg, pvor_cen
     38  REAL,    PARAMETER :: sg0=0.75  !--- Start level for PV=cte search loop
     39  INTEGER, PARAMETER :: nadj=3    !--- Adjacent levs nb for thresholds detection
     40  REAL,    PARAMETER :: w(5)=[0.1,0.25,0.3,0.25,0.1] !--- Vertical smoothing
     41  INTEGER, SAVE :: k0
     42  LOGICAL, SAVE :: first=.TRUE.
     43!$OMP THREADPRIVATE(k0,first)
    3944!-------------------------------------------------------------------------------
    4045  sub='dyn_tropopause'
    … …  
    4247  CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
    4348  CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
    44   CALL assert(SIZE(presnivs,1)==klev, TRIM(sub)//" presnivs klev")
    4549  CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
    4650  CALL assert(SHAPE(pplay)==[klon,klev  ],TRIM(sub)//" pplay shape")
    … …  
    4852
    4953  !--- DEFAULT THRESHOLDS
    50   tpo0=380.;   IF(PRESENT(tpot0)) tpo0=tpot0  !--- In kelvins
    51   pvo0=2.0;    IF(PRESENT(pvor0)) pvo0=pvor0  !--- In PVU
    52   pmn0= 8000.; IF(PRESENT(pmin0)) pmn0=pmin0  !--- In pascals
    53   pmx0=50000.; IF(PRESENT(pmax0)) pmx0=pmax0  !--- In pascals
    54   kmin=klev-locate(presnivs(klev:1:-1),pmx0)+1
    55   kmax=klev-locate(presnivs(klev:1:-1),pmn0)+1
     54  th0=380.; IF(PRESENT(thet0)) th0=thet0   !--- In kelvins
     55  pv0=  2.; IF(PRESENT(pvor0)) pv0=pvor0   !--- In PVU
     56  IF(first) THEN
     57    DO k0=1,klev; IF(aps(k0)/preff+bps(k0)<sg0) EXIT; END DO; first=.FALSE.
     58  END IF
    5659
    57   !--- POTENTIAL TEMPERATURE AT CELLS CENTERS
    58   DO k= 1, klev
    59     DO i = 1, klon
    60       t_pot(i,k) = t(i,k)*(preff/pplay(i,k))**RKAPPA
     60  !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES
     61  DO i = 1,klon
     62    tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA
     63    tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA
     64    DO k=2,klev
     65      al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
     66      tpot_cen(i,k) =  t(i,k)                        *(preff/pplay(i,k))**RKAPPA
     67      tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA
     68      !--- FORCE QUANTITIES TO BE GROWING
     69      IF(tpot_edg(i,k)<tpot_edg(i,k-1)) tpot_edg(i,k)=tpot_edg(i,k-1)+1.E-5
     70      IF(tpot_cen(i,k)<tpot_cen(i,k-1)) tpot_cen(i,k)=tpot_cen(i,k-1)+1.E-5
    6171    END DO
    62   END DO
    63 
    64   !--- TEMPERATURE AT CELLS INTERFACES (except in top layer)
    65   t_edg(:,1) = ts(:)
    66   DO k= 2, klev
    67     DO i = 1, klon
    68       al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
    69       t_edg(i,k) = t(i,k) + al*(t(i,k-1)-t(i,k))
    70     END DO
     72    !--- VERTICAL SMOOTHING
     73    tpot_cen(i,:)=smooth(tpot_cen(i,:),w)
     74    tpot_edg(i,:)=smooth(tpot_edg(i,:),w)
    7175  END DO
    7276
    7377  !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
    74   DO k= 1, klev-1
    75     DO i = 1, klon
    76       IF(paprs(i,k)==paprs(i,k+1)) THEN; pvort(i,k)=HUGE(1.); CYCLE; END IF
    77       pvort(i,k) = -1.E6 * RG * 2.*ROMEGA*ABS(SIN(latitude_deg(i)*RPI/180.))   &
    78                          * ( t_edg(i,k+1)*(preff/paprs(i,k+1) )**RKAPPA        &
    79                            - t_edg(i,k  )*(preff/paprs(i,k  ) )**RKAPPA)       &
    80                          / ( paprs(i,k+1)  -     paprs(i,k  ) )
     78  DO i = 1, klon
     79    DO k= 1, klev-1
     80      pvor_cen(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))&
     81                   * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k))
    8182    END DO
     83    !--- VERTICAL SMOOTHING
     84    pvor_cen(i,1:klev-1)=smooth(pvor_cen(i,1:klev-1),w)
    8285  END DO
    8386
    84   !--- LOCATE TROPOPAUSE
    85   ptrop(:)=0.
     87  CALL writefield_phy('PV_ou'   ,pvor_cen(:,1:klev-1),klev-1)
     88  CALL writefield_phy('theta_ou',tpot_cen,klev)
     89  CALL writefield_phy('Pcen_ou' ,pplay   ,klev)
     90
     91  !--- LOCATE TROPOPAUSE: LOWEST POINT BETWEEN THETA=380K AND PV=2PVU SURFACES.
    8692  DO i = 1, klon
    87     !--- Dynamical tropopause
    88     it=kmax; DO WHILE(pvort(i,it)>pvo0.AND.it>=kmin); it=it-1; END DO
    89     IF(pvort(i,it+1)/=pvort(i,it).AND.ALL([kmax,kmin-1]/=it) &
    90       .AND.ALL([pvort(i,it),pvort(i,it+1)]/=HUGE(1.))) THEN
    91       al = (pvort(i,it+1)-pvo0)/(pvort(i,it+1)-pvort(i,it))
    92       ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
    93     END IF
    94     !--- Potential temperature iso-surface
    95     it = kmin-1+locate(t_pot(i,kmin:kmax),tpo0)
    96     IF(t_pot(i,it+1)/=t_pot(i,it).AND.ALL([kmin-1,kmax]/=it)) THEN
    97       al = (t_pot(i,it+1)-tpo0)/(t_pot(i,it+1)-t_pot(i,it))
    98       ptrop(i) = MAX(ptrop(i),pplay(i,it)**al * pplay(i,it+1)**(1.-al))
    99     END IF
     93    !--- UPPER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM TOP
     94    DO kt=klev-1,1,-1; IF(ALL(ABS(pvor_cen(i,kt-nadj:kt))<=pv0)) EXIT; END DO
     95    !--- LOWER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM BOTTOM
     96    DO kb=k0,klev-1;   IF(ALL(ABS(pvor_cen(i,kb:kb+nadj))> pv0)) EXIT; END DO; kb=kb-1
     97    !--- ISO-THETA POINT: THETA=380K       STARTING FROM TOP
     98    DO kp=klev-1,1,-1; IF(ALL(ABS(tpot_cen(i,kp-nadj:kp))<=th0)) EXIT; END DO
     99    !--- CHOOSE BETWEEN LOWER AND UPPER TROPOPAUSE
     100    IF(2*COUNT(ABS(pvor_cen(i,kb:kt))>pv0)>kt-kb+1) kt=kb
     101    !--- PV-DEFINED TROPOPAUSE
     102    al = (ABS(pvor_cen(i,kt+1))-pv0)/ABS(pvor_cen(i,kt+1)-pvor_cen(i,kt))
     103    dyn_tropopause(i) = pplay(i,kt+1)*(pplay(i,kt)/pplay(i,kt+1))**al
     104    !--- THETA=380K IN THE TROPICAL REGION
     105    al = (tpot_cen(i,kp+1)-th0)/(tpot_cen(i,kp+1)-tpot_cen(i,kp))
     106    dyn_tropopause(i) = MAX( pplay(i,kp+1)*(pplay(i,kp)/pplay(i,kp+1))**al,    &
     107                            dyn_tropopause(i) )
     108  END DO
     109  CALL writefield_phy('PreTr_ou',dyn_tropopause,1)
     110
     111END FUNCTION dyn_tropopause
     112
     113
     114!-------------------------------------------------------------------------------
     115!
     116FUNCTION smooth(v,w)
     117!
     118!-------------------------------------------------------------------------------
     119! Arguments:
     120  REAL, INTENT(IN)         :: v(:), w(:)
     121  REAL, DIMENSION(SIZE(v)) :: smooth
     122!-------------------------------------------------------------------------------
     123! Local variables:
     124  INTEGER :: nv, nw, k, kb, ke, lb, le
     125!-------------------------------------------------------------------------------
     126  nv=SIZE(v); nw=(SIZE(w)-1)/2
     127  DO k=1,nv
     128    kb=MAX(k-nw,1 ); lb=MAX(2+nw   -k,1)
     129    ke=MIN(k+nw,nv); le=MIN(1+nw+nv-k,1+2*nw)
     130    smooth(k)=SUM(v(kb:ke)*w(lb:le))/SUM(w(lb:le))
    100131  END DO
    101132
    102 END SUBROUTINE dyn_tropopause
     133END FUNCTION smooth
     134!
     135!-------------------------------------------------------------------------------
    103136
    104137END MODULE tropopause_m
Note: See TracChangeset for help on using the changeset viewer.