Index: LMDZ6/trunk/libf/misc/regr_conserv_m.F90 =================================================================== --- LMDZ6/trunk/libf/misc/regr_conserv_m.F90 (revision 3068) +++ LMDZ6/trunk/libf/misc/regr_conserv_m.F90 (revision 3069) @@ -84,5 +84,5 @@ IF(a==xt(it )) co=co+xt(it )-xs(is ) IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1) - vt(it) = vt(it)+idt*(b-a)*(vs(is)+co*slope(is)) + vt(it) = vt(it)+idt*(b-a)*(vs(is)+co*slope(is)/2.) ELSE vt(it) = vt(it)+idt*(b-a)* vs(is) @@ -142,6 +142,6 @@ IF(a==xt(it )) co=co+xt(it )-xs(is ) IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1) - IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)*(vs(is,:)+co*slope(is,:)) - IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)*(vs(:,is)+co*slope(:,is)) + IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)*(vs(is,:)+co*slope(is,:)/2.) + IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)*(vs(:,is)+co*slope(:,is)/2.) ELSE IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)* vs(is,:) @@ -202,7 +202,7 @@ IF(a==xt(it )) co=co+xt(it )-xs(is ) IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1) - IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)*(vs(is,:,:)+co*slope(is,:,:)) - IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)*(vs(:,is,:)+co*slope(:,is,:)) - IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)*(vs(:,:,is)+co*slope(:,:,is)) + IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)*(vs(is,:,:)+co*slope(is,:,:)/2.) + IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)*(vs(:,is,:)+co*slope(:,is,:)/2.) + IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)*(vs(:,:,is)+co*slope(:,:,is)/2.) ELSE IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)* vs(is,:,:) @@ -264,8 +264,8 @@ IF(a==xt(it )) co=co+xt(it )-xs(is ) IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1) - IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)*(vs(is,:,:,:)+co*slope(is,:,:,:)) - IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)*(vs(:,is,:,:)+co*slope(:,is,:,:)) - IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)*(vs(:,:,is,:)+co*slope(:,:,is,:)) - IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)*(vs(:,:,:,is)+co*slope(:,:,:,is)) + IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)*(vs(is,:,:,:)+co*slope(is,:,:,:)/2.) + IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)*(vs(:,is,:,:)+co*slope(:,is,:,:)/2.) + IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)*(vs(:,:,is,:)+co*slope(:,:,is,:)/2.) + IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)*(vs(:,:,:,is)+co*slope(:,:,:,is)/2.) ELSE IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)* vs(is,:,:,:) @@ -328,9 +328,9 @@ IF(a==xt(it )) co=co+xt(it )-xs(is ) IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1) - IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)*(vs(is,:,:,:,:)+co*slope(is,:,:,:,:)) - IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)*(vs(:,is,:,:,:)+co*slope(:,is,:,:,:)) - IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)*(vs(:,:,is,:,:)+co*slope(:,:,is,:,:)) - IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)*(vs(:,:,:,is,:)+co*slope(:,:,:,is,:)) - IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)*(vs(:,:,:,:,is)+co*slope(:,:,:,:,is)) + IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)*(vs(is,:,:,:,:)+co*slope(is,:,:,:,:)/2.) + IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)*(vs(:,is,:,:,:)+co*slope(:,is,:,:,:)/2.) + IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)*(vs(:,:,is,:,:)+co*slope(:,:,is,:,:)/2.) + IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)*(vs(:,:,:,is,:)+co*slope(:,:,:,is,:)/2.) + IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)*(vs(:,:,:,:,is)+co*slope(:,:,:,:,is)/2.) ELSE IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)* vs(is,:,:,:,:) Index: LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90 (revision 3068) +++ LMDZ6/trunk/libf/phylmd/regr_pr_time_av_m.F90 (revision 3069) @@ -5,4 +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,9 +87,14 @@ 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,16 +123,17 @@ !------------------------------------------------------------------------------- ! 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,4 +165,7 @@ 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,21 +264,30 @@ 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_ouSigT_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,6 +301,7 @@ 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,5 +311,5 @@ *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,5 +327,30 @@ ! 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,4 +593,5 @@ 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