Changeset 3069 for LMDZ6/trunk/libf/phylmd

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.

File:

• LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90 (modified) (9 diffs)

LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90

@@ -5 +5 @@
   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
 
@@ -86 +87 @@
   LOGICAL, SAVE :: lfirst=.TRUE.                   !--- First call flag
 !$OMP THREADPRIVATE(lfirst)
-  REAL,    PARAMETER :: pTropUp=5.E+3 !--- Value < tropopause pressure (Pa)
+  REAL,    PARAMETER :: pTropUp=9.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
+  LOGICAL, PARAMETER :: debug=.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
@@ -117 +123 @@
 !-------------------------------------------------------------------------------
 ! Arguments:
-  INTEGER, INTENT(IN)  :: fID           !--- NetCDF file ID
+  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)  :: julien              !--- Days since Jan 1st
+  !--- File & LMDZ (resp. decreasing & increasing order) interfaces pressure, Pa
+  REAL,    INTENT(IN)  :: pint_in(:)          !--- in:  file          (nlev_in+1)
+  REAL,    INTENT(IN)  :: pint_ou(:,:)        !--- out: LMDZ     (klon,nbp_lev+1)
+  REAL,    INTENT(OUT) :: v3(:,:,:)           !--- Regridded fld (klon,llm,n_var)
   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)
+  REAL,    INTENT(IN), OPTIONAL :: lon_in(:)  !--- File longitudes vector (klon)
+  REAL,    INTENT(IN), OPTIONAL :: lat_in(:)  !--- File latitudes  vector (klon)
+  REAL,    INTENT(IN), OPTIONAL :: pcen_in(:) !--- File mid-layers pres (nlev_in)
+  REAL,    INTENT(IN), OPTIONAL :: ptrop_ou(:)!--- LMDZ tropopause pres   (klon)
 !-------------------------------------------------------------------------------
 ! Local variables:
@@ -158 +165 @@
   REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
   LOGICAL :: ll
+!--- For debug
+  REAL, DIMENSION(klon)           :: ptrop_in, ptrop_ef!, dzStrain, dzStrain0
+! REAL, DIMENSION(klon,nbp_lev+1) :: pstrn_ou, phii
 !-------------------------------------------------------------------------------
   sub="regr_pr_time_av"
@@ -254 +264 @@
       Sig_ou(:) = [pintou (1:nbp_lev)/ps2(i),1.0] !--- increasing values
 
-      !--- INPUT (FILE) AND OUTPUT (LMDZ) SIGMA LEVELS AT TROPOPAUSE
+      !--- 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.
       SigT_in = get_SigTrop(i,itrp)
+      SigT_in=MIN(SigT_in,ChemPTrMax/ps2(i))      !--- too low  value filtered
+      SigT_in=MAX(SigT_in,ChemPTrMin/ps2(i))      !--- too high value filtered
+
+      !--- 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.
       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
+      IF(SigT_ou<SigT_in/rho) SigT_ou=SigT_in/rho !--- too low  value w/r input
+      IF(SigT_ou>SigT_in*rho) SigT_ou=SigT_in*rho !--- too high value w/r input
+      SigT_ou=MIN(SigT_ou,DynPTrMax/ps2(i))       !--- too low  value filtered
+      SigT_ou=MAX(SigT_ou,DynPTrMin/ps2(i))       !--- too high value filtered
       ptropou(i)=SigT_ou*ps2(i)
 
-      !--- STRETCHING EXPONENT INCREMENT FOR SIMPLE POWER LAW
+      !--- FULLY STRETCHED LAYER BOUNDS (alpha power law fully applied)
       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
+      !--- PARTIALLY STRETCHED LAYER BOUNDS (fraction of alpha applied)
+      ! The used formulae ensure the first derivative stays positive.
       beta = LOG(Sig_top)/LOG(Sig_bot)
       Sig_bo0 = Sig_bot ; IF(alpha<0.) Sig_bo0 = Sig_bot**(1/beta)
@@ -282 +301 @@
       ito0 = locate(Sig_ou(:),Sig_to0)
 
-      !--- FUNCTION FOR STRETCHING LOCALISATION
-      !    0 < Sig_to0 < Sig_top <= Sig_bo0 < Sig_bot < 1
+      !--- STRETCHING POWER LAW FUNCTION:
+      !    0    in [0,Sig_to0] and [Sig_bo0,1] ;    1    in [Sig_bot,Sig_top]
+      !    0->1 in [Sig_to0,Sig_top]                1->0 in [Sig_bot,Sig_bo0]
       phi(:)=0.
       phi(itop+1:ibot) =  1.
@@ -291 +311 @@
                             *LOG(Sig_bot)/LOG(Sig_bot/Sig_bo0)
 
-      !--- LOCAL STRAINED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
+      !--- LOCALY STRECHED OUTPUT (LMDZ) PRESSURE PROFILES (INCREASING ORDER)
       pstr_ou(:) = pintou(:) * Sig_ou(:)**(alpha*phi(:))
 
@@ -307 +327 @@
 !     IF(ll) CALL abort_physic(sub, 'Inconsistent O3 values in stratosphere', 1)
 
+      IF(debug) THEN
+!        dzStrain0(i)=7.*LOG(Sig_bo0/Sig_to0)
+!        dzStrain (i)=7.*LOG(Sig_bot/Sig_top)
+        ptrop_in(i) = SigT_in*ps2(i)
+!        Pstrn_ou(i,:)= pstr_ou
+!        phii(i,:)=phi(:)
+        ptrop_ef(i)=chem_tropopause(i, itrp, locate(pintou(:),PTropUp),        &
+                           pintou(:), v3(:,nbp_lev+1:1:-1,1),o3trop=o3t0)
+      END IF
     END DO
+    IF(debug) THEN
+!      CALL writefield_phy('PreSt_ou' ,Pstrn_ou,nbp_lev+1) !--- Strained pressure
+!      CALL writefield_phy('dzStrain' ,dzStrain ,1)   !--- Fully & total strained
+!      CALL writefield_phy('dzStrain0',dzStrain0,1)   !--- regions thickness
+!      CALL writefield_phy('phi',phii,nbp_lev+1)      !--- Localisation function
+      !--- Tropopauses pressures:
+      CALL writefield_phy('PreTr_in',ptrop_in,1)     !--- Input and effective
+      CALL writefield_phy('PreTr_ef',ptrop_ef,1)     !    chemical tropopauses
+    END IF
+  END IF
+  IF(debug) THEN
+    IF(lAdjTro) &
+    CALL writefield_phy('PreTr_ou',ptropou,1)        !--- LMDz dyn tropopause
+    CALL writefield_phy('Ozone_in',v2(:,:,1),nlev_in)!--- Raw input O3 field
+    CALL writefield_phy('Ozone_ou',v3(:,:,1),nbp_lev)!--- Output ozone field
+    CALL writefield_phy('Pint_ou' ,pint_ou,1+nbp_lev)!--- LMDZ unstreched Press
   END IF
 
@@ -548 +593 @@
   lmax=.FALSE.; IF(PRESENT(vmax)) lmax=COUNT(o3col>vmax)/=0
   check_ozone=lmin.OR.lmax; IF(.NOT.check_ozone) RETURN
+  IF(prt_level<1) RETURN
 
   !--- SOME TOO LOW VALUES FOUND