Changeset 3069


Ignore:
Timestamp:
Nov 13, 2017, 6:11:11 PM (7 years ago)
Author:
dcugnet
Message:

Fix the ozone vertical interpolation routines. Main features:

  • lower initial pressure for chemical tropopause detection routine,

to avoid ozone holes, in 1998 in particular.

  • min and max tropopauses pressure ensured for both field (chemical)

and lmdz (dynamical) tropopauses pressures for more safety.

  • fix in the second order term of interpolation.
Location:
LMDZ6/trunk/libf
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/misc/regr_conserv_m.F90

    r2788 r3069  
    8484    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    8585    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    86     vt(it) = vt(it)+idt*(b-a)*(vs(is)+co*slope(is))
     86    vt(it) = vt(it)+idt*(b-a)*(vs(is)+co*slope(is)/2.)
    8787  ELSE
    8888    vt(it) = vt(it)+idt*(b-a)* vs(is)
    … …  
    142142    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    143143    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    144     IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)*(vs(is,:)+co*slope(is,:))
    145     IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)*(vs(:,is)+co*slope(:,is))
     144    IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)*(vs(is,:)+co*slope(is,:)/2.)
     145    IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)*(vs(:,is)+co*slope(:,is)/2.)
    146146  ELSE
    147147    IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)* vs(is,:)
    … …  
    202202    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    203203    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    204     IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)*(vs(is,:,:)+co*slope(is,:,:))
    205     IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)*(vs(:,is,:)+co*slope(:,is,:))
    206     IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)*(vs(:,:,is)+co*slope(:,:,is))
     204    IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)*(vs(is,:,:)+co*slope(is,:,:)/2.)
     205    IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)*(vs(:,is,:)+co*slope(:,is,:)/2.)
     206    IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)*(vs(:,:,is)+co*slope(:,:,is)/2.)
    207207  ELSE
    208208    IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)* vs(is,:,:)
    … …  
    264264    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    265265    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    266     IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)*(vs(is,:,:,:)+co*slope(is,:,:,:))
    267     IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)*(vs(:,is,:,:)+co*slope(:,is,:,:))
    268     IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)*(vs(:,:,is,:)+co*slope(:,:,is,:))
    269     IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)*(vs(:,:,:,is)+co*slope(:,:,:,is))
     266    IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)*(vs(is,:,:,:)+co*slope(is,:,:,:)/2.)
     267    IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)*(vs(:,is,:,:)+co*slope(:,is,:,:)/2.)
     268    IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)*(vs(:,:,is,:)+co*slope(:,:,is,:)/2.)
     269    IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)*(vs(:,:,:,is)+co*slope(:,:,:,is)/2.)
    270270  ELSE
    271271    IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)* vs(is,:,:,:)
    … …  
    328328    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    329329    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    330     IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)*(vs(is,:,:,:,:)+co*slope(is,:,:,:,:))
    331     IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)*(vs(:,is,:,:,:)+co*slope(:,is,:,:,:))
    332     IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)*(vs(:,:,is,:,:)+co*slope(:,:,is,:,:))
    333     IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)*(vs(:,:,:,is,:)+co*slope(:,:,:,is,:))
    334     IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)*(vs(:,:,:,:,is)+co*slope(:,:,:,:,is))
     330    IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)*(vs(is,:,:,:,:)+co*slope(is,:,:,:,:)/2.)
     331    IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)*(vs(:,is,:,:,:)+co*slope(:,is,:,:,:)/2.)
     332    IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)*(vs(:,:,is,:,:)+co*slope(:,:,is,:,:)/2.)
     333    IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)*(vs(:,:,:,is,:)+co*slope(:,:,:,is,:)/2.)
     334    IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)*(vs(:,:,:,:,is)+co*slope(:,:,:,:,is)/2.)
    335335  ELSE
    336336    IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)* vs(is,:,:,:,:)

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

    r2981 r3069  
    55  USE mod_phys_lmdz_transfert_para, ONLY: bcast
    66  USE mod_phys_lmdz_para, ONLY: mpi_rank, omp_rank
     7  USE print_control_mod,  ONLY: prt_level
    78  IMPLICIT NONE
    89
    … …  
    8687  LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
    8788!$OMP THREADPRIVATE(lfirst)
    88   REAL,    PARAMETER :: pTropUp=5.E+3 !--- Value < tropopause pressure (Pa)
     89  REAL,    PARAMETER :: pTropUp=9.E+3 !--- Value  < tropopause pressure (Pa)
    8990  REAL,    PARAMETER :: gamm = 0.4    !--- Relative thickness of transitions
    9091  REAL,    PARAMETER :: rho  = 1.3    !--- Max tropopauses sigma ratio
    9192  REAL,    PARAMETER :: o3t0 = 1.E-7  !--- Nominal O3 vmr at tropopause
    9293  LOGICAL, PARAMETER :: lo3tp=.FALSE. !--- Use parametrized O3 vmr at tropopause
     94  LOGICAL, PARAMETER :: debug=.FALSE. !--- Force writefield_phy multiple outputs
     95  REAL, PARAMETER :: ChemPTrMin=9.E+3 !--- Thresholds for minimum and maximum
     96  REAL, PARAMETER :: ChemPTrMax=4.E+4 !    chemical  tropopause pressure (Pa).
     97  REAL, PARAMETER :: DynPTrMin =8.E+3 !--- Thresholds for minimum and maximum
     98  REAL, PARAMETER :: DynPTrMax =4.E+4 !    dynamical tropopause pressure (Pa).
    9399
    94100CONTAINS
    … …  
    117123!-------------------------------------------------------------------------------
    118124! Arguments:
    119   INTEGER, INTENT(IN)  :: fID           !--- NetCDF file ID
     125  INTEGER, INTENT(IN)  :: fID                 !--- NetCDF file ID
    120126  CHARACTER(LEN=13), INTENT(IN) :: nam(:)     !--- NetCDF variables names
    121   REAL,    INTENT(IN)  :: julien        !--- Days since Jan 1st
    122   REAL,    INTENT(IN)  :: pint_in(:)    !--- Interfaces file Pres, Pa, ascending
    123   REAL,    INTENT(IN)  :: pint_ou(:,:)  !--- Interfaces LMDZ pressure, Pa (g,nbp_lev+1)
    124   REAL,    INTENT(OUT) :: v3(:,:,:)     !--- Regridded field (klon,llm,SIZE(nam))
     127  REAL,    INTENT(IN)  :: julien              !--- Days since Jan 1st
     128  !--- File & LMDZ (resp. decreasing & increasing order) interfaces pressure, Pa
     129  REAL,    INTENT(IN)  :: pint_in(:)          !--- in:  file          (nlev_in+1)
     130  REAL,    INTENT(IN)  :: pint_ou(:,:)        !--- out: LMDZ     (klon,nbp_lev+1)
     131  REAL,    INTENT(OUT) :: v3(:,:,:)           !--- Regridded fld (klon,llm,n_var)
    125132  REAL,    INTENT(IN), OPTIONAL :: time_in(:) !--- Records times, in days
    126133                                              !    since Jan 1 of current year
    127   REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- Input/output longitudes vector
    128   REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- Input/output latitudes vector
    129   REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- Mid-layers file pressure
    130   REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- Output tropopause pres (klon)
     134  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- File longitudes vector (klon)
     135  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- File latitudes  vector (klon)
     136  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- File mid-layers pres (nlev_in)
     137  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- LMDZ tropopause pres  (klon)
    131138!-------------------------------------------------------------------------------
    132139! Local variables:
    … …  
    158165  REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
    159166  LOGICAL :: ll
     167!--- For debug
     168  REAL, DIMENSION(klon)           :: ptrop_in, ptrop_ef!, dzStrain, dzStrain0
     169! REAL, DIMENSION(klon,nbp_lev+1) :: pstrn_ou, phii
    160170!-------------------------------------------------------------------------------
    161171  sub="regr_pr_time_av"
    … …  
    254264      Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values
    255265
    256       !--- INPUT (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE
     266      !--- INPUT (FILE) SIGMA LEVEL AT TROPOPAUSE ; extreme values are filtered
     267      ! to keep tropopause pressure realistic ; high values are usually due to
     268      ! ozone hole fooling the crude chemical tropopause detection algorithm.
    257269      SigT_in = get_SigTrop(i,itrp)
     270      SigT_in=MIN(SigT_in,ChemPTrMax/ps2(i))      !--- too low  value filtered
     271      SigT_in=MAX(SigT_in,ChemPTrMin/ps2(i))      !--- too high value filtered
     272
     273      !--- OUTPUT (LMDZ) SIGMA LEVEL AT TROPOPAUSE ; too high variations of the
     274      ! dynamical tropopause (especially in filaments) are conterbalanced with a
     275      ! filter ensuring it stays within a certain distance around input (file)
     276      ! tropopause, hence avoiding avoid a too thick stretched region ; a final
     277      ! extra-safety filter keeps the tropopause pressure value realistic.
    258278      SigT_ou = ptrop_ou(i)/ps2(i)
    259 
    260       !--- AVOID THE FILAMENTS WHICH WOULD NEED A VERY THICK STRETCHED REGION
    261       IF(SigT_ou>SigT_in*rho) SigT_ou = SigT_in*rho
    262       IF(SigT_ou<SigT_in/rho) SigT_ou = SigT_in/rho
     279      IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho !--- too low  value w/r input
     280      IF(SigT_ou>SigT_in*rho) SigT_ou=SigT_in*rho !--- too high value w/r input
     281      SigT_ou=MIN(SigT_ou,DynPTrMax/ps2(i))       !--- too low  value filtered
     282      SigT_ou=MAX(SigT_ou,DynPTrMin/ps2(i))       !--- too high value filtered
    263283      ptropou(i)=SigT_ou*ps2(i)
    264284
    265       !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW
     285      !--- FULLY STRETCHED LAYER BOUNDS (alpha power law fully applied)
    266286      alpha = LOG(SigT_in/SigT_ou)/LOG(SigT_ou)
    267 
    268       !--- FULLY STRETCHED LAYER BOUNDS (FILE AND MODEL TROPOPAUSES)
    269287      Sig_bot = MAX(SigT_in,SigT_ou) ; ibot = locate(Sig_ou(:),Sig_bot)
    270288      Sig_top = MIN(SigT_in,SigT_ou) ; itop = locate(Sig_ou(:),Sig_top)
    271289
    272       !--- PARTIALLY STRETCHED LAYER BOUNDS, ENSURING >0 DERIVATIVE
     290      !--- PARTIALLY STRETCHED LAYER BOUNDS (fraction of alpha applied)
     291      ! The used formulae ensure the first derivative stays positive.
    273292      beta = LOG(Sig_top)/LOG(Sig_bot)
    274293      Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
    … …  
    282301      ito0 = locate(Sig_ou(:),Sig_to0)
    283302
    284       !--- FUNCTION FOR STRETCHING LOCALISATION
    285       !    0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1
     303      !--- STRETCHING POWER LAW FUNCTION:
     304      !    0    in [0,Sig_to0] and [Sig_bo0,1] ;    1    in [Sig_bot,Sig_top]
     305      !    0->1 in [Sig_to0,Sig_top]                1->0 in [Sig_bot,Sig_bo0]
    286306      phi(:)=0.
    287307      phi(itop+1:ibot) =  1.
    … …  
    291311                            *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)
    292312
    293       !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
     313      !--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
    294314      pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
    295315
    … …  
    307327!     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
    308328
     329      IF(debug) THEN
     330!        dzStrain0(i)=7.*LOG(Sig_bo0/Sig_to0)
     331!        dzStrain (i)=7.*LOG(Sig_bot/Sig_top)
     332        ptrop_in(i) = SigT_in*ps2(i)
     333!        Pstrn_ou(i,:)= pstr_ou
     334!        phii(i,:)=phi(:)
     335        ptrop_ef(i)=chem_tropopause(i, itrp, locate(pintou(:),PTropUp),        &
     336                           pintou(:), v3(:,nbp_lev+1:1:-1,1),o3trop=o3t0)
     337      END IF
    309338    END DO
     339    IF(debug) THEN
     340!      CALL writefield_phy('PreSt_ou' ,Pstrn_ou,nbp_lev+1) !--- Strained pressure
     341!      CALL writefield_phy('dzStrain' ,dzStrain ,1)   !--- Fully & total strained
     342!      CALL writefield_phy('dzStrain0',dzStrain0,1)   !--- regions thickness
     343!      CALL writefield_phy('phi',phii,nbp_lev+1)      !--- Localisation function
     344      !--- Tropopauses pressures:
     345      CALL writefield_phy('PreTr_in',ptrop_in,1)     !--- Input and effective
     346      CALL writefield_phy('PreTr_ef',ptrop_ef,1)     !    chemical tropopauses
     347    END IF
     348  END IF
     349  IF(debug) THEN
     350    IF(lAdjTro) &
     351    CALL writefield_phy('PreTr_ou',ptropou,1)        !--- LMDz dyn tropopause
     352    CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
     353    CALL writefield_phy('Ozone_ou',v3(:,:,1),nbp_lev)!--- Output ozone field
     354    CALL writefield_phy('Pint_ou' ,pint_ou,1+nbp_lev)!--- LMDZ unstreched Press
    310355  END IF
    311356
    … …  
    548593  lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
    549594  check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN
     595  IF(prt_level<1) RETURN
    550596
    551597  !--- SOME TOO LOW VALUES FOUND
Note: See TracChangeset for help on using the changeset viewer.