Context Navigation

lmdz_lscp.F90 @ 5220

Last change on this file since 5220 was 5194, checked in by abarral, 2 months ago
Merge r5189, r5191 from trunk
File size: 63.1 KB

Rev	Line
[4664]	1	MODULE lmdz_lscp
[3999]	2
[5111]	3	IMPLICIT NONE
[3999]	4
	5	CONTAINS
	6
[5111]	7	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	8	SUBROUTINE lscp(klon, klev, dtime, missing_val, &
	9	paprs, pplay, temp, qt, qice_save, ptconv, ratqs, &
	10	d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, &
	11	pfraclr, pfracld, &
	12	cldfraliq, sigma2_icefracturb, mean_icefracturb, &
	13	radocond, radicefrac, rain, snow, &
	14	frac_impa, frac_nucl, beta, &
	15	prfl, psfl, rhcl, qta, fraca, &
	16	tv, pspsk, tla, thl, iflag_cld_th, &
	17	iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, &
	18	distcltop, temp_cltop, tke, tke_dissip, &
	19	qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
	20	Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
	21	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv, &
	22	qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, &
	23	dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim, &
	24	dqsmelt, dqsfreez)
[3999]	25
[5111]	26	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	27	! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
	28	! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
	29	!--------------------------------------------------------------------------------
	30	! Date: 01/2021
	31	!--------------------------------------------------------------------------------
	32	! Aim: Large Scale Clouds and Precipitation (LSCP)
[5099]	33
[5111]	34	! This code is a new version of the fisrtilp.F90 routine, which is itself a
	35	! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
	36	! and 'ilp' (il pleut, L. Li)
[5099]	37
[5111]	38	! Compared to the original fisrtilp code, lscp
	39	! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
	40	! -> consider always precipitation thermalisation (fl_cor_ebil>0)
	41	! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
	42	! -> option "oldbug" by JYG has been removed
	43	! -> iflag_t_glace >0 always
	44	! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
	45	! -> rectangular distribution from L. Li no longer available
	46	! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
	47	!--------------------------------------------------------------------------------
	48	! References:
[5099]	49
[5111]	50	! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
	51	! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
	52	! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
	53	! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
	54	! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
	55	! - Touzze-Peifert Ludo, PhD thesis, p117-124
	56	! -------------------------------------------------------------------------------
	57	! Code structure:
[5099]	58
[5111]	59	! P0> Thermalization of the precipitation coming from the overlying layer
	60	! P1> Evaporation of the precipitation (falling from the k+1 level)
	61	! P2> Cloud formation (at the k level)
	62	! P2.A.1> With the PDFs, calculation of cloud properties using the inital
	63	! values of T and Q
	64	! P2.A.2> Coupling between condensed water and temperature
	65	! P2.A.3> Calculation of final quantities associated with cloud formation
	66	! P2.B> Release of Latent heat after cloud formation
	67	! P3> Autoconversion to precipitation (k-level)
	68	! P4> Wet scavenging
	69	!------------------------------------------------------------------------------
	70	! Some preliminary comments (JBM) :
[5099]	71
[5111]	72	! The cloud water that the radiation scheme sees is not the same that the cloud
	73	! water used in the physics and the dynamics
[5099]	74
[5111]	75	! During the autoconversion to precipitation (P3 step), radocond (cloud water used
	76	! by the radiation scheme) is calculated as an average of the water that remains
	77	! in the cloud during the precipitation and not the water remaining at the end
	78	! of the time step. The latter is used in the rest of the physics and advected
	79	! by the dynamics.
[5099]	80
[5111]	81	! In summary:
[5099]	82
[5111]	83	! Radiation:
	84	! xflwc(newmicro)+xfiwc(newmicro) =
	85	! radocond=lwcon(nc)+iwcon(nc)
[5099]	86
[5111]	87	! Notetheless, be aware of:
[5099]	88
[5111]	89	! radocond .NE. ocond(nc)
	90	! i.e.:
	91	! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
	92	! but oliq+(ocond-oliq) .EQ. ocond
	93	! (which is not trivial)
	94	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4654]	95
[5111]	96	! USE de modules contenant des fonctions.
	97	USE lmdz_cloudth, ONLY: cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
	98	USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat
	99	USE lmdz_lscp_tools, ONLY: icefrac_lscp, icefrac_lscp_turb
	100	USE lmdz_lscp_tools, ONLY: fallice_velocity, distance_to_cloud_top
	101	USE ice_sursat_mod, ONLY: ice_sursat
	102	USE lmdz_lscp_poprecip, ONLY: poprecip_precld, poprecip_postcld
[3999]	103
[5111]	104	! Use du module lmdz_lscp_ini contenant les constantes
	105	USE lmdz_lscp_ini, ONLY: prt_level, lunout
	106	USE lmdz_lscp_ini, ONLY: seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
	107	USE lmdz_lscp_ini, ONLY: ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
	108	USE lmdz_lscp_ini, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
	109	USE lmdz_lscp_ini, ONLY: coef_eva, coef_sub, cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
	110	USE lmdz_lscp_ini, ONLY: iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
	111	USE lmdz_lscp_ini, ONLY: iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
	112	USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
	113	USE lmdz_lscp_ini, ONLY: ok_poprecip
	114	USE lmdz_lscp_ini, ONLY: iflag_icefrac
	115	USE lmdz_abort_physic, ONLY: abort_physic
	116	IMPLICIT NONE
[3999]	117
[5111]	118	!===============================================================================
	119	! VARIABLES DECLARATION
	120	!===============================================================================
[3999]	121
[5111]	122	! INPUT VARIABLES:
	123	!-----------------
[3999]	124
[5111]	125	INTEGER, INTENT(IN) :: klon, klev ! number of horizontal grid points and vertical levels
	126	REAL, INTENT(IN) :: dtime ! time step [s]
	127	REAL, INTENT(IN) :: missing_val ! missing value for output
[4059]	128
[5111]	129	REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs ! inter-layer pressure [Pa]
	130	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa]
	131	REAL, DIMENSION(klon, klev), INTENT(IN) :: temp ! temperature (K)
	132	REAL, DIMENSION(klon, klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg]
	133	REAL, DIMENSION(klon, klev), INTENT(IN) :: qice_save ! ice specific from previous time step [kg/kg]
	134	INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds
	135	INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
[5116]	136	! CR: if iflag_ice_thermo=2, ONLY convection is active
[5111]	137	LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated
[4059]	138
[5111]	139	LOGICAL, DIMENSION(klon, klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
[3999]	140
[5111]	141	!Inputs associated with thermal plumes
[4226]	142
[5111]	143	REAL, DIMENSION(klon, klev), INTENT(IN) :: tv ! virtual potential temperature [K]
	144	REAL, DIMENSION(klon, klev), INTENT(IN) :: qta ! specific humidity within thermals [kg/kg]
	145	REAL, DIMENSION(klon, klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-]
	146	REAL, DIMENSION(klon, klev), INTENT(IN) :: pspsk ! exner potential (p/100000)**(R/cp)
	147	REAL, DIMENSION(klon, klev), INTENT(IN) :: tla ! liquid temperature within thermals [K]
	148	REAL, DIMENSION(klon, klev), INTENT(IN) :: tke !--turbulent kinetic energy [m2/s2]
	149	REAL, DIMENSION(klon, klev), INTENT(IN) :: tke_dissip !--TKE dissipation [m2/s3]
[3999]	150
[5111]	151	! INPUT/OUTPUT variables
	152	!------------------------
[4059]	153
[5111]	154	REAL, DIMENSION(klon, klev), INTENT(INOUT) :: thl ! liquid potential temperature [K]
	155	REAL, DIMENSION(klon, klev), INTENT(INOUT) :: ratqs ! function of pressure that sets the large-scale
[4380]	156
[3999]	157
[5111]	158	! Input sursaturation en glace
	159	REAL, DIMENSION(klon, klev), INTENT(INOUT) :: rneb_seri ! fraction nuageuse en memoire
[4686]	160
[5111]	161	! OUTPUT variables
	162	!-----------------
[3999]	163
[5111]	164	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! temperature increment [K]
	165	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg]
	166	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg]
	167	REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg]
	168	REAL, DIMENSION(klon, klev), INTENT(OUT) :: rneb ! cloud fraction [-]
	169	REAL, DIMENSION(klon, klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-]
	170	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-]
	171	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-]
	172	REAL, DIMENSION(klon, klev), INTENT(OUT) :: cldfraliq ! liquid fraction of cloud [-]
	173	REAL, DIMENSION(klon, klev), INTENT(OUT) :: sigma2_icefracturb ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
	174	REAL, DIMENSION(klon, klev), INTENT(OUT) :: mean_icefracturb ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
	175	REAL, DIMENSION(klon, klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg]
	176	REAL, DIMENSION(klon, klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme
	177	REAL, DIMENSION(klon, klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-]
	178	REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2]
	179	REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2]
	180	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg]
	181	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg]
	182	REAL, DIMENSION(klon, klev + 1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2]
	183	REAL, DIMENSION(klon, klev + 1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2]
	184	REAL, DIMENSION(klon, klev), INTENT(OUT) :: distcltop ! distance to cloud top [m]
	185	REAL, DIMENSION(klon, klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K]
	186	REAL, DIMENSION(klon, klev), INTENT(OUT) :: beta ! conversion rate of condensed water
[4686]	187
[5111]	188	! fraction of aerosol scavenging through impaction and nucleation (for on-line)
[4686]	189
[5111]	190	REAL, DIMENSION(klon, klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-]
	191	REAL, DIMENSION(klon, klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-]
[4803]	192
[5111]	193	! for supersaturation and contrails parameterisation
[3999]	194
[5111]	195	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qclr ! specific total water content in clear sky region [kg/kg]
	196	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qcld ! specific total water content in cloudy region [kg/kg]
	197	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qss ! specific total water content in supersat region [kg/kg]
	198	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg]
	199	REAL, DIMENSION(klon, klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-]
	200	REAL, DIMENSION(klon, klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-]
	201	REAL, DIMENSION(klon, klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-]
	202	REAL, DIMENSION(klon, klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K]
	203	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg]
	204	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q)
	205	REAL, DIMENSION(klon, klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails
	206	REAL, DIMENSION(klon, klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails
	207	REAL, DIMENSION(klon, klev), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals
	208	REAL, DIMENSION(klon, klev), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment
	209	REAL, DIMENSION(klon, klev), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals
	210	REAL, DIMENSION(klon, klev), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment
[3999]	211
[5111]	212	! for POPRECIP
[3999]	213
[5111]	214	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg]
	215	REAL, DIMENSION(klon, klev), INTENT(OUT) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg]
	216	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqreva !--rain tendendy due to evaporation [kg/kg/s]
	217	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
	218	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
	219	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
	220	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
	221	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
	222	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqsrim !--snow tendency due to riming [kg/kg/s]
	223	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqsmelt !--snow tendency due to melting [kg/kg/s]
	224	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqrmelt !--rain tendency due to melting [kg/kg/s]
	225	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqsfreez !--snow tendency due to freezing [kg/kg/s]
	226	REAL, DIMENSION(klon, klev), INTENT(OUT) :: dqrfreez !--rain tendency due to freezing [kg/kg/s]
[3999]	227
	228
[5111]	229	! LOCAL VARIABLES:
	230	!----------------
	231	REAL, DIMENSION(klon) :: qsl, qsi ! saturation threshold at current vertical level
	232	REAL :: zct, zcl, zexpo
	233	REAL, DIMENSION(klon, klev) :: ctot
	234	REAL, DIMENSION(klon, klev) :: ctot_vol
	235	REAL, DIMENSION(klon) :: zqs, zdqs
	236	REAL :: zdelta, zcor, zcvm5
	237	REAL, DIMENSION(klon) :: zdqsdT_raw
	238	REAL, DIMENSION(klon) :: gammasat, dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T
	239	REAL, DIMENSION(klon) :: Tbef, qlbef, DT ! temperature, humidity and temp. variation during lognormal iteration
	240	REAL :: num, denom
	241	REAL :: cste
	242	REAL, DIMENSION(klon) :: zpdf_sig, zpdf_k, zpdf_delta ! lognormal parameters
	243	REAL, DIMENSION(klon) :: Zpdf_a, zpdf_b, zpdf_e1, zpdf_e2 ! lognormal intermediate variables
	244	REAL :: erf
	245	REAL, DIMENSION(klon) :: zfice_th
	246	REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
	247	REAL, DIMENSION(klon) :: zrfl, zrfln
	248	REAL :: zqev, zqevt
	249	REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
	250	REAL :: zqev0, zqevi, zqevti
	251	REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
	252	REAL, DIMENSION(klon) :: zoliql, zoliqi
	253	REAL, DIMENSION(klon) :: zt
	254	REAL, DIMENSION(klon, klev) :: zrho
	255	REAL, DIMENSION(klon) :: zdz, iwc
	256	REAL :: zchau, zfroi
	257	REAL, DIMENSION(klon) :: zfice, zneb, znebprecip
	258	REAL :: zmelt, zrain, zsnow, zprecip
	259	REAL, DIMENSION(klon) :: dzfice
	260	REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
	261	REAL :: zsolid
	262	REAL, DIMENSION(klon) :: qtot, qzero
	263	REAL, DIMENSION(klon) :: dqsl, dqsi
	264	REAL :: smallestreal
	265	! Variables for Bergeron process
	266	REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
	267	REAL, DIMENSION(klon) :: zqpreci, zqprecl
	268	! Variables precipitation energy conservation
	269	REAL, DIMENSION(klon) :: zmqc
	270	REAL :: zalpha_tr
	271	REAL :: zfrac_lessi
	272	REAL, DIMENSION(klon) :: zprec_cond
	273	REAL :: zmair
	274	REAL :: zcpair, zcpeau
	275	REAL, DIMENSION(klon) :: zlh_solid
	276	REAL, DIMENSION(klon) :: ztupnew
	277	REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
	278	REAL :: zm_solid ! for liquid -> solid conversion
	279	REAL, DIMENSION(klon) :: zrflclr, zrflcld
	280	REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
	281	REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
	282	REAL, DIMENSION(klon) :: ziflclr, ziflcld
	283	REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
	284	REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
	285	REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
	286	REAL, DIMENSION(klon, klev) :: velo
	287	REAL :: vr, ffallv
	288	REAL :: qlmpc, qimpc, rnebmpc
	289	REAL, DIMENSION(klon, klev) :: radocondi, radocondl
	290	REAL :: effective_zneb
	291	REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
	292	REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl ! for icefrac_lscp_turb
[3999]	293
[5111]	294	INTEGER i, k, n, kk, iter
	295	INTEGER, DIMENSION(klon) :: n_i
	296	INTEGER ncoreczq
	297	INTEGER, DIMENSION(klon, klev) :: mpc_bl_points
	298	LOGICAL iftop
[4654]	299
[5111]	300	LOGICAL, DIMENSION(klon) :: lognormale
	301	LOGICAL, DIMENSION(klon) :: keepgoing
[3999]	302
[5111]	303	CHARACTER (len = 20) :: modname = 'lscp'
	304	CHARACTER (len = 80) :: abort_message
[3999]	305
	306
[5111]	307	!===============================================================================
	308	! INITIALISATION
	309	!===============================================================================
[3999]	310
[5111]	311	! Few initial checks
[3999]	312
[5111]	313	IF (iflag_fisrtilp_qsat < 0) THEN
	314	abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
	315	CALL abort_physic(modname, abort_message, 1)
	316	ENDIF
[4380]	317
[5111]	318	! Few initialisations
[5007]	319
[5111]	320	znebprecip(:) = 0.0
	321	ctot_vol(1:klon, 1:klev) = 0.0
	322	rneblsvol(1:klon, 1:klev) = 0.0
	323	smallestreal = 1.e-9
	324	znebprecipclr(:) = 0.0
	325	znebprecipcld(:) = 0.0
	326	mpc_bl_points(:, :) = 0
[3999]	327
[5111]	328	IF (prt_level>9) WRITE(lunout, *) 'NUAGES4 A. JAM'
[4803]	329
[5111]	330	! AA for 'safety' reasons
	331	zalpha_tr = 0.
	332	zfrac_lessi = 0.
	333	beta(:, :) = 0.
[4803]	334
[5111]	335	! Initialisation of variables:
[4392]	336
[5111]	337	prfl(:, :) = 0.0
	338	psfl(:, :) = 0.0
	339	d_t(:, :) = 0.0
	340	d_q(:, :) = 0.0
	341	d_ql(:, :) = 0.0
	342	d_qi(:, :) = 0.0
	343	rneb(:, :) = 0.0
	344	pfraclr(:, :) = 0.0
	345	pfracld(:, :) = 0.0
	346	cldfraliq(:, :) = 0.
	347	sigma2_icefracturb(:, :) = 0.
	348	mean_icefracturb(:, :) = 0.
	349	radocond(:, :) = 0.0
	350	radicefrac(:, :) = 0.0
	351	frac_nucl(:, :) = 1.0
	352	frac_impa(:, :) = 1.0
	353	rain(:) = 0.0
	354	snow(:) = 0.0
	355	zoliq(:) = 0.0
	356	zfice(:) = 0.0
	357	dzfice(:) = 0.0
	358	zfice_turb(:) = 0.0
	359	dzfice_turb(:) = 0.0
	360	zqprecl(:) = 0.0
	361	zqpreci(:) = 0.0
	362	zrfl(:) = 0.0
	363	zifl(:) = 0.0
	364	ziflprev(:) = 0.0
	365	zneb(:) = seuil_neb
	366	zrflclr(:) = 0.0
	367	ziflclr(:) = 0.0
	368	zrflcld(:) = 0.0
	369	ziflcld(:) = 0.0
	370	tot_zneb(:) = 0.0
	371	tot_znebn(:) = 0.0
	372	d_tot_zneb(:) = 0.0
	373	qzero(:) = 0.0
	374	zdistcltop(:) = 0.0
	375	ztemp_cltop(:) = 0.0
	376	ztupnew(:) = 0.0
	377	!--ice supersaturation
	378	gamma_ss(:, :) = 1.0
	379	qss(:, :) = 0.0
	380	rnebss(:, :) = 0.0
	381	Tcontr(:, :) = missing_val
	382	qcontr(:, :) = missing_val
	383	qcontr2(:, :) = missing_val
	384	fcontrN(:, :) = 0.0
	385	fcontrP(:, :) = 0.0
	386	distcltop(:, :) = 0.
	387	temp_cltop(:, :) = 0.
[4392]	388
[5111]	389	!-- poprecip
	390	qraindiag(:, :) = 0.
	391	qsnowdiag(:, :) = 0.
	392	dqreva(:, :) = 0.
	393	dqrauto(:, :) = 0.
	394	dqrmelt(:, :) = 0.
	395	dqrfreez(:, :) = 0.
	396	dqrcol(:, :) = 0.
	397	dqssub(:, :) = 0.
	398	dqsauto(:, :) = 0.
	399	dqsrim(:, :) = 0.
	400	dqsagg(:, :) = 0.
	401	dqsfreez(:, :) = 0.
	402	dqsmelt(:, :) = 0.
	403	zqupnew(:) = 0.
	404	zqvapclr(:) = 0.
[3999]	405
	406
	407
[5111]	408	!c_iso: variable initialisation for iso
[4803]	409
[5111]	410
	411	!===============================================================================
	412	! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
	413	!===============================================================================
	414
	415	ncoreczq = 0
	416
	417	DO k = klev, 1, -1
	418
	419	IF (k<=klev - 1) THEN
	420	iftop = .FALSE.
	421	ELSE
	422	iftop = .TRUE.
	423	ENDIF
	424
	425	! Initialisation temperature and specific humidity
	426	! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
	427	! at the end of the klon loop, a temperature incremtent d_t due to all processes
	428	! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
	429	! d_t = temperature tendency due to lscp
	430	! The temperature of the overlying layer is updated here because needed for thermalization
	431	DO i = 1, klon
	432	zt(i) = temp(i, k)
	433	zq(i) = qt(i, k)
[5117]	434	IF (.NOT. iftop) THEN
[5111]	435	ztupnew(i) = temp(i, k + 1) + d_t(i, k + 1)
	436	zqupnew(i) = qt(i, k + 1) + d_q(i, k + 1) + d_ql(i, k + 1) + d_qi(i, k + 1)
	437	!--zqs(i) is the saturation specific humidity in the layer above
	438	zqvapclr(i) = MAX(0., qt(i, k + 1) + d_q(i, k + 1) - rneb(i, k + 1) * zqs(i))
[4803]	439	ENDIF
[4392]	440	!c_iso init of iso
[5111]	441	ENDDO
[4803]	442
[5111]	443	!================================================================
	444	! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
	445	IF (ok_poprecip) THEN
[4803]	446
[5111]	447	CALL poprecip_precld(klon, dtime, iftop, paprs(:, k), paprs(:, k + 1), pplay(:, k), &
	448	zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
	449	zqvapclr, zqupnew, &
	450	zrfl, zrflclr, zrflcld, &
	451	zifl, ziflclr, ziflcld, &
	452	dqreva(:, k), dqssub(:, k) &
	453	)
[4803]	454
[5111]	455	!================================================================
	456	ELSE
[3999]	457
[5111]	458	! --------------------------------------------------------------------
	459	! P1> Thermalization of precipitation falling from the overlying layer
	460	! --------------------------------------------------------------------
	461	! Computes air temperature variation due to enthalpy transported by
	462	! precipitation. Precipitation is then thermalized with the air in the
	463	! layer.
	464	! The precipitation should remain thermalized throughout the different
	465	! thermodynamical transformations.
	466	! The corresponding water mass should
	467	! be added when calculating the layer's enthalpy change with
	468	! temperature
	469	! See lmdzpedia page todoan
	470	! todoan: check consistency with ice phase
	471	! todoan: understand why several steps
	472	! ---------------------------------------------------------------------
	473
	474	IF (iftop) THEN
	475
	476	DO i = 1, klon
[4803]	477	zmqc(i) = 0.
[5111]	478	ENDDO
[4803]	479
[5111]	480	ELSE
[4803]	481
[5111]	482	DO i = 1, klon
	483
	484	zmair = (paprs(i, k) - paprs(i, k + 1)) / RG
[3999]	485	! no condensed water so cp=cp(vapor+dry air)
	486	! RVTMP2=rcpv/rcpd-1
[5111]	487	zcpair = RCPD * (1.0 + RVTMP2 * zq(i))
	488	zcpeau = RCPD * RVTMP2
[4226]	489
[5111]	490	! zmqc: precipitation mass that has to be thermalized with
[3999]	491	! layer's air so that precipitation at the ground has the
	492	! same temperature as the lowermost layer
[5111]	493	zmqc(i) = (zrfl(i) + zifl(i)) * dtime / zmair
[3999]	494	! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
[5111]	495	zt(i) = (ztupnew(i) * zmqc(i) * zcpeau + zcpair * zt(i)) &
	496	/ (zcpair + zmqc(i) * zcpeau)
[4226]	497
[5111]	498	ENDDO
[3999]	499
[5111]	500	ENDIF
[3999]	501
[5111]	502	! --------------------------------------------------------------------
	503	! P2> Precipitation evaporation/sublimation/melting
	504	! --------------------------------------------------------------------
	505	! A part of the precipitation coming from above is evaporated/sublimated/melted.
	506	! --------------------------------------------------------------------
[3999]	507
[5111]	508	IF (iflag_evap_prec>=1) THEN
[4226]	509
[5111]	510	! Calculation of saturation specific humidity
	511	! depending on temperature:
	512	CALL calc_qsat_ecmwf(klon, zt(:), qzero(:), pplay(:, k), RTT, 0, .FALSE., zqs(:), zdqs(:))
	513	! wrt liquid water
	514	CALL calc_qsat_ecmwf(klon, zt(:), qzero(:), pplay(:, k), RTT, 1, .FALSE., qsl(:), dqsl(:))
	515	! wrt ice
	516	CALL calc_qsat_ecmwf(klon, zt(:), qzero(:), pplay(:, k), RTT, 2, .FALSE., qsi(:), dqsi(:))
	517
	518	DO i = 1, klon
	519
[3999]	520	! if precipitation
[5111]	521	IF (zrfl(i) + zifl(i)>0.) THEN
[4226]	522
[5111]	523	! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
	524	! c_iso: likely important to distinguish cs from neb isotope precipitation
[4392]	525
[5111]	526	IF (iflag_evap_prec>=4) THEN
	527	zrfl(i) = zrflclr(i)
	528	zifl(i) = ziflclr(i)
	529	ENDIF
[3999]	530
[5111]	531	IF (iflag_evap_prec==1) THEN
	532	znebprecip(i) = zneb(i)
	533	ELSE
	534	znebprecip(i) = MAX(zneb(i), znebprecip(i))
	535	ENDIF
[3999]	536
[5111]	537	IF (iflag_evap_prec>4) THEN
[4563]	538	! Max evaporation not to saturate the clear sky precip fraction
	539	! i.e. the fraction where evaporation occurs
[5111]	540	zqev0 = MAX(0.0, (zqs(i) - zq(i)) * znebprecipclr(i))
	541	ELSEIF (iflag_evap_prec == 4) THEN
[4226]	542	! Max evaporation not to saturate the whole mesh
[4563]	543	! Pay attention -> lead to unrealistic and excessive evaporation
[5111]	544	zqev0 = MAX(0.0, zqs(i) - zq(i))
	545	ELSE
[4226]	546	! Max evap not to saturate the fraction below the cloud
[5111]	547	zqev0 = MAX(0.0, (zqs(i) - zq(i)) * znebprecip(i))
	548	ENDIF
[3999]	549
[5111]	550	! Evaporation of liquid precipitation coming from above
	551	! dP/dz=beta(1-q/qsat)sqrt(P)
	552	! formula from Sundquist 1988, Klemp & Wilhemson 1978
	553	! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
[3999]	554
[5111]	555	IF (iflag_evap_prec==3) THEN
	556	zqevt = znebprecip(i) * coef_eva * (1.0 - zq(i) / qsl(i)) &
	557	* SQRT(zrfl(i) / max(1.e-4, znebprecip(i))) &
	558	* (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG
	559	ELSE IF (iflag_evap_prec>=4) THEN
	560	zqevt = znebprecipclr(i) * coef_eva * (1.0 - zq(i) / qsl(i)) &
	561	* SQRT(zrfl(i) / max(1.e-8, znebprecipclr(i))) &
	562	* (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG
	563	ELSE
	564	zqevt = 1. * coef_eva * (1.0 - zq(i) / qsl(i)) * SQRT(zrfl(i)) &
	565	* (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG
	566	ENDIF
[3999]	567
[5111]	568	zqevt = MAX(0.0, MIN(zqevt, zrfl(i))) &
	569	* RG * dtime / (paprs(i, k) - paprs(i, k + 1))
[3999]	570
[5111]	571	! sublimation of the solid precipitation coming from above
	572	IF (iflag_evap_prec==3) THEN
	573	zqevti = znebprecip(i) * coef_sub * (1.0 - zq(i) / qsi(i)) &
	574	* SQRT(zifl(i) / max(1.e-4, znebprecip(i))) &
	575	* (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG
	576	ELSE IF (iflag_evap_prec>=4) THEN
	577	zqevti = znebprecipclr(i) * coef_sub * (1.0 - zq(i) / qsi(i)) &
	578	* SQRT(zifl(i) / max(1.e-8, znebprecipclr(i))) &
	579	* (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG
	580	ELSE
	581	zqevti = 1. * coef_sub * (1.0 - zq(i) / qsi(i)) * SQRT(zifl(i)) &
	582	* (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * RD / RG
	583	ENDIF
[4226]	584
[5111]	585	zqevti = MAX(0.0, MIN(zqevti, zifl(i))) &
	586	* RG * dtime / (paprs(i, k) - paprs(i, k + 1))
[3999]	587
[5111]	588	! A. JAM
	589	! Evaporation limit: we ensure that the layer's fraction below
	590	! the cloud or the whole mesh (depending on iflag_evap_prec)
	591	! does not reach saturation. In this case, we
	592	! redistribute zqev0 conserving the ratio liquid/ice
[4226]	593
[5111]	594	IF (zqevt + zqevti>zqev0) THEN
	595	zqev = zqev0 * zqevt / (zqevt + zqevti)
	596	zqevi = zqev0 * zqevti / (zqevt + zqevti)
	597	ELSE
	598	zqev = zqevt
	599	zqevi = zqevti
	600	ENDIF
[3999]	601
	602
[5111]	603	! New solid and liquid precipitation fluxes after evap and sublimation
	604	zrfln(i) = Max(0., zrfl(i) - zqev * (paprs(i, k) - paprs(i, k + 1)) &
	605	/ RG / dtime)
	606	zifln(i) = Max(0., zifl(i) - zqevi * (paprs(i, k) - paprs(i, k + 1)) &
	607	/ RG / dtime)
[3999]	608
[4392]	609
[5111]	610	! vapor, temperature, precip fluxes update
	611	! vapor is updated after evaporation/sublimation (it is increased)
	612	zq(i) = zq(i) - (zrfln(i) + zifln(i) - zrfl(i) - zifl(i)) &
	613	* (RG / (paprs(i, k) - paprs(i, k + 1))) * dtime
	614	! zmqc is the total condensed water in the precip flux (it is decreased)
	615	zmqc(i) = zmqc(i) + (zrfln(i) + zifln(i) - zrfl(i) - zifl(i)) &
	616	* (RG / (paprs(i, k) - paprs(i, k + 1))) * dtime
	617	! air and precip temperature (i.e., gridbox temperature)
	618	! is updated due to latent heat cooling
	619	zt(i) = zt(i) + (zrfln(i) - zrfl(i)) &
	620	* (RG / (paprs(i, k) - paprs(i, k + 1))) * dtime &
	621	* RLVTT / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i))) &
	622	+ (zifln(i) - zifl(i)) &
	623	* (RG / (paprs(i, k) - paprs(i, k + 1))) * dtime &
	624	* RLSTT / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i)))
[3999]	625
[5111]	626	! New values of liquid and solid precipitation
	627	zrfl(i) = zrfln(i)
	628	zifl(i) = zifln(i)
[3999]	629
[5111]	630	! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
	631	! due to evap + sublim
[4392]	632
[5111]	633	IF (iflag_evap_prec>=4) THEN
	634	zrflclr(i) = zrfl(i)
	635	ziflclr(i) = zifl(i)
	636	IF(zrflclr(i) + ziflclr(i)<=0) THEN
	637	znebprecipclr(i) = 0.0
[4226]	638	ENDIF
[5111]	639	zrfl(i) = zrflclr(i) + zrflcld(i)
	640	zifl(i) = ziflclr(i) + ziflcld(i)
	641	ENDIF
[3999]	642
[5111]	643	! c_iso duplicate for isotopes or loop on isotopes
[3999]	644
[5111]	645	! Melting:
	646	zmelt = ((zt(i) - RTT) / (ztfondue - RTT)) ! JYG
	647	! precip fraction that is melted
	648	zmelt = MIN(MAX(zmelt, 0.), 1.)
[4392]	649
[5111]	650	! update of rainfall and snowfall due to melting
	651	IF (iflag_evap_prec>=4) THEN
	652	zrflclr(i) = zrflclr(i) + zmelt * ziflclr(i)
	653	zrflcld(i) = zrflcld(i) + zmelt * ziflcld(i)
	654	zrfl(i) = zrflclr(i) + zrflcld(i)
	655	ELSE
	656	zrfl(i) = zrfl(i) + zmelt * zifl(i)
	657	ENDIF
[4412]	658
	659
[5111]	660	! c_iso: melting of isotopic precipi with zmelt (no fractionation)
[3999]	661
[5111]	662	! Latent heat of melting because of precipitation melting
	663	! NB: the air + precip temperature is simultaneously updated
	664	zt(i) = zt(i) - zifl(i) * zmelt * (RG * dtime) / (paprs(i, k) - paprs(i, k + 1)) &
	665	* RLMLT / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i)))
[3999]	666
[5111]	667	IF (iflag_evap_prec>=4) THEN
	668	ziflclr(i) = ziflclr(i) * (1. - zmelt)
	669	ziflcld(i) = ziflcld(i) * (1. - zmelt)
	670	zifl(i) = ziflclr(i) + ziflcld(i)
	671	ELSE
	672	zifl(i) = zifl(i) * (1. - zmelt)
	673	ENDIF
	674
[3999]	675	ELSE
[5111]	676	! if no precip, we reinitialize the cloud fraction used for the precip to 0
	677	znebprecip(i) = 0.
[3999]	678
	679	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	680
[5111]	681	ENDDO ! loop on klon
[4226]	682
[5111]	683	ENDIF ! (iflag_evap_prec>=1)
[3999]	684
[5111]	685	ENDIF ! (ok_poprecip)
[4803]	686
[5111]	687	! --------------------------------------------------------------------
	688	! End precip evaporation
	689	! --------------------------------------------------------------------
[4226]	690
[5111]	691	! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
	692	!-------------------------------------------------------
[3999]	693
[5111]	694	qtot(:) = zq(:) + zmqc(:)
	695	CALL calc_qsat_ecmwf(klon, zt(:), qtot(:), pplay(:, k), RTT, 0, .FALSE., zqs(:), zdqs(:))
	696	DO i = 1, klon
	697	zdelta = MAX(0., SIGN(1., RTT - zt(i)))
	698	zdqsdT_raw(i) = zdqs(i) * RCPD * (1.0 + RVTMP2 * zq(i)) / (RLVTT * (1. - zdelta) + RLSTT * zdelta)
	699	IF (zq(i) < 1.e-15) THEN
	700	ncoreczq = ncoreczq + 1
	701	zq(i) = 1.e-15
	702	ENDIF
	703	! c_iso: do something similar for isotopes
[5099]	704
[5111]	705	ENDDO
[3999]	706
[5111]	707	! --------------------------------------------------------------------
	708	! P2> Cloud formation
	709	!---------------------------------------------------------------------
[3999]	710
[5111]	711	! Unlike fisrtilp, we always assume a 'fractional cloud' approach
	712	! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
	713	! is prescribed and depends on large scale variables and boundary layer
	714	! properties)
	715	! The decrease in condensed part due tu latent heating is taken into
	716	! account
	717	! -------------------------------------------------------------------
[3999]	718
[5111]	719	! P2.1> With the PDFs (log-normal, bigaussian)
	720	! cloud properties calculation with the initial values of t and q
	721	! ----------------------------------------------------------------
[3999]	722
[5111]	723	! initialise gammasat and qincloud_mpc
	724	gammasat(:) = 1.
	725	qincloud_mpc(:) = 0.
[3999]	726
[5111]	727	IF (iflag_cld_th>=5) THEN
	728	! Cloud cover and content in meshes affected by shallow convection,
	729	! are retrieved from a bi-gaussian distribution of the saturation deficit
	730	! following Jam et al. 2013
[3999]	731
[5111]	732	IF (iflag_cloudth_vert<=2) THEN
	733	! Old version of Arnaud Jam
[4651]	734
[5111]	735	CALL cloudth(klon, klev, k, tv, &
	736	zq, qta, fraca, &
	737	qcloud, ctot, pspsk, paprs, pplay, tla, thl, &
	738	ratqs, zqs, temp, &
	739	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[3999]	740
[5111]	741	ELSEIF (iflag_cloudth_vert>=3 .AND. iflag_cloudth_vert<=5) THEN
	742	! Default version of Arnaud Jam
[3999]	743
[5111]	744	CALL cloudth_v3(klon, klev, k, tv, &
	745	zq, qta, fraca, &
	746	qcloud, ctot, ctot_vol, pspsk, paprs, pplay, tla, thl, &
	747	ratqs, zqs, temp, &
	748	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[3999]	749
[5111]	750	ELSEIF (iflag_cloudth_vert==6) THEN
	751	! Jean Jouhaud's version, with specific separation between surface and volume
	752	! cloud fraction Decembre 2018
[3999]	753
[5111]	754	CALL cloudth_v6(klon, klev, k, tv, &
	755	zq, qta, fraca, &
	756	qcloud, ctot, ctot_vol, pspsk, paprs, pplay, tla, thl, &
	757	ratqs, zqs, temp, &
	758	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[4910]	759
[5111]	760	ELSEIF (iflag_cloudth_vert == 7) THEN
	761	! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
	762	! for boundary-layer mixed phase clouds
	763	CALL cloudth_mpc(klon, klev, k, mpc_bl_points, zt, zq, qta(:, k), fraca(:, k), &
	764	pspsk(:, k), paprs(:, k + 1), paprs(:, k), pplay(:, k), tla(:, k), &
	765	ratqs(:, k), qcloud, qincloud_mpc, zfice_th, ctot(:, k), ctot_vol(:, k), &
	766	cloudth_sth(:, k), cloudth_senv(:, k), cloudth_sigmath(:, k), cloudth_sigmaenv(:, k))
[3999]	767
[5111]	768	ENDIF
[3999]	769
[5111]	770	DO i = 1, klon
	771	rneb(i, k) = ctot(i, k)
	772	rneblsvol(i, k) = ctot_vol(i, k)
	773	zqn(i) = qcloud(i)
	774	ENDDO
[3999]	775
[5111]	776	ENDIF
[3999]	777
[5111]	778	IF (iflag_cld_th <= 4) THEN
[3999]	779
[5111]	780	! lognormal
	781	lognormale(:) = .TRUE.
[3999]	782
[5111]	783	ELSEIF (iflag_cld_th >= 6) THEN
[3999]	784
[5111]	785	! lognormal distribution when no thermals
	786	lognormale(:) = fraca(:, k) < min_frac_th_cld
[3999]	787
[5111]	788	ELSE
	789	! When iflag_cld_th=5, we always assume
	790	! bi-gaussian distribution
	791	lognormale(:) = .FALSE.
[4226]	792
[5111]	793	ENDIF
[4226]	794
[5111]	795	DT(:) = 0.
	796	n_i(:) = 0
	797	Tbef(:) = zt(:)
	798	qlbef(:) = 0.
[3999]	799
[5111]	800	! Treatment of non-boundary layer clouds (lognormale)
	801	! condensation with qsat(T) variation (adaptation)
	802	! Iterative resolution to converge towards qsat
	803	! with update of temperature, ice fraction and qsat at
	804	! each iteration
[4226]	805
[5111]	806	! todoan -> sensitivity to iflag_fisrtilp_qsat
	807	DO iter = 1, iflag_fisrtilp_qsat + 1
[4226]	808
[5111]	809	DO i = 1, klon
[3999]	810
[5111]	811	! keepgoing = .TRUE. while convergence is not satisfied
	812	keepgoing(i) = ABS(DT(i))>DDT0
[3999]	813
[5111]	814	IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN
[3999]	815
[5111]	816	! if not convergence:
[4226]	817
[5111]	818	! P2.2.1> cloud fraction and condensed water mass calculation
	819	! Calculated variables:
	820	! rneb : cloud fraction
	821	! zqn : total water within the cloud
	822	! zcond: mean condensed water within the mesh
	823	! rhcl: clear-sky relative humidity
	824	!---------------------------------------------------------------
[3999]	825
[5111]	826	! new temperature that only serves in the iteration process:
	827	Tbef(i) = Tbef(i) + DT(i)
[4072]	828
[5111]	829	! Rneb, qzn and zcond for lognormal PDFs
	830	qtot(i) = zq(i) + zmqc(i)
[4072]	831
[5111]	832	ENDIF
[5007]	833
[5111]	834	ENDDO
[4072]	835
[5111]	836	! Calculation of saturation specific humidity and ice fraction
	837	CALL calc_qsat_ecmwf(klon, Tbef(:), qtot(:), pplay(:, k), RTT, 0, .FALSE., zqs(:), zdqs(:))
	838	CALL calc_gammasat(klon, Tbef(:), qtot(:), pplay(:, k), gammasat(:), dgammasatdt(:))
	839	! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
	840	zdqs(:) = gammasat(:) * zdqs(:) + zqs(:) * dgammasatdt(:)
	841	! cloud phase determination
	842	IF (iflag_t_glace>=4) THEN
	843	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
	844	CALL distance_to_cloud_top(klon, klev, k, temp, pplay, paprs, rneb, zdistcltop, ztemp_cltop)
	845	ENDIF
[4072]	846
[5111]	847	CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:), ztemp_cltop(:), zfice(:), dzfice(:))
[4226]	848
[5111]	849	DO i = 1, klon !todoan : check if loop in i is needed
[4059]	850
[5111]	851	IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN
[4059]	852
[5111]	853	zpdf_sig(i) = ratqs(i, k) * zq(i)
	854	zpdf_k(i) = -sqrt(log(1. + (zpdf_sig(i) / zq(i))**2))
	855	zpdf_delta(i) = log(zq(i) / (gammasat(i) * zqs(i)))
	856	zpdf_a(i) = zpdf_delta(i) / (zpdf_k(i) * sqrt(2.))
	857	zpdf_b(i) = zpdf_k(i) / (2. * sqrt(2.))
	858	zpdf_e1(i) = zpdf_a(i) - zpdf_b(i)
	859	zpdf_e1(i) = sign(min(ABS(zpdf_e1(i)), 5.), zpdf_e1(i))
	860	zpdf_e1(i) = 1. - erf(zpdf_e1(i))
	861	zpdf_e2(i) = zpdf_a(i) + zpdf_b(i)
	862	zpdf_e2(i) = sign(min(ABS(zpdf_e2(i)), 5.), zpdf_e2(i))
	863	zpdf_e2(i) = 1. - erf(zpdf_e2(i))
[4226]	864
[5111]	865	IF ((.NOT.ok_ice_sursat).OR.(Tbef(i)>t_glace_min)) THEN
[4226]	866
[5111]	867	IF (zpdf_e1(i)<1.e-10) THEN
	868	rneb(i, k) = 0.
	869	zqn(i) = gammasat(i) * zqs(i)
	870	ELSE
	871	rneb(i, k) = 0.5 * zpdf_e1(i)
	872	zqn(i) = zq(i) * zpdf_e2(i) / zpdf_e1(i)
	873	ENDIF
[3999]	874
[5111]	875	rnebss(i, k) = 0.0 !--added by OB (needed because of the convergence loop on time)
	876	fcontrN(i, k) = 0.0 !--idem
	877	fcontrP(i, k) = 0.0 !--idem
	878	qss(i, k) = 0.0 !--idem
[4059]	879
[5111]	880	ELSE ! in case of ice supersaturation by Audran
[4059]	881
[5111]	882	!------------------------------------
	883	! ICE SUPERSATURATION
	884	!------------------------------------
[3999]	885
[5111]	886	CALL ice_sursat(pplay(i, k), paprs(i, k) - paprs(i, k + 1), dtime, i, k, temp(i, k), zq(i), &
	887	gamma_ss(i, k), zqs(i), Tbef(i), rneb_seri(i, k), ratqs(i, k), &
	888	rneb(i, k), zqn(i), rnebss(i, k), qss(i, k), &
	889	Tcontr(i, k), qcontr(i, k), qcontr2(i, k), fcontrN(i, k), fcontrP(i, k))
[3999]	890
[5117]	891	ENDIF ! ((flag_ice_sursat.EQ.0).OR.(Tbef(i).gt.t_glace_min))
[3999]	892
[5111]	893	! If vertical heterogeneity, change fraction by volume as well
	894	IF (iflag_cloudth_vert>=3) THEN
	895	ctot_vol(i, k) = rneb(i, k)
	896	rneblsvol(i, k) = ctot_vol(i, k)
	897	ENDIF
[3999]	898
	899
[5111]	900	! P2.2.2> Approximative calculation of temperature variation DT
	901	! due to condensation.
	902	! Calculated variables:
	903	! dT : temperature change due to condensation
	904	!---------------------------------------------------------------
[3999]	905
[5111]	906	IF (zfice(i)<1) THEN
	907	cste = RLVTT
	908	ELSE
	909	cste = RLSTT
	910	ENDIF
[3999]	911
[5111]	912	! LEA_R : check formule
	913	qlbef(i) = max(0., zqn(i) - zqs(i))
	914	num = -Tbef(i) + zt(i) + rneb(i, k) * ((1 - zfice(i)) * RLVTT &
	915	+ zfice(i) * RLSTT) / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i))) * qlbef(i)
	916	denom = 1. + rneb(i, k) * ((1 - zfice(i)) * RLVTT + zfice(i) * RLSTT) / cste * zdqs(i) &
	917	- (RLSTT - RLVTT) / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i))) * rneb(i, k) &
	918	* qlbef(i) * dzfice(i)
	919	! here we update a provisory temperature variable that only serves in the iteration
	920	! process
	921	DT(i) = num / denom
	922	n_i(i) = n_i(i) + 1
[3999]	923
[5111]	924	ENDIF ! end keepgoing
[3999]	925
[5111]	926	ENDDO ! end loop on i
[3999]	927
[5111]	928	ENDDO ! iter=1,iflag_fisrtilp_qsat+1
[4562]	929
[5111]	930	! P2.3> Final quantities calculation
	931	! Calculated variables:
	932	! rneb : cloud fraction
	933	! zcond: mean condensed water in the mesh
	934	! zqn : mean water vapor in the mesh
	935	! zfice: ice fraction in clouds
	936	! zt : temperature
	937	! rhcl : clear-sky relative humidity
	938	! ----------------------------------------------------------------
[5007]	939
	940
[5111]	941	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
	942	IF (iflag_t_glace>=4) THEN
	943	CALL distance_to_cloud_top(klon, klev, k, temp, pplay, paprs, rneb, zdistcltop, ztemp_cltop)
	944	distcltop(:, k) = zdistcltop(:)
	945	temp_cltop(:, k) = ztemp_cltop(:)
	946	ENDIF
[5007]	947
[5111]	948	! Partition function depending on temperature
[5194]	949	CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
[3999]	950
[5111]	951	! Partition function depending on tke for non shallow-convective clouds
	952	IF (iflag_icefrac >= 1) THEN
[3999]	953
[5111]	954	CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:, k), paprs(:, k), paprs(:, k + 1), qice_save(:, k), ziflcld, zqn, &
	955	rneb(:, k), tke(:, k), tke_dissip(:, k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:, k), sigma2_icefracturb(:, k), mean_icefracturb(:, k))
	956
	957	ENDIF
	958
	959	! Water vapor update, Phase determination and subsequent latent heat exchange
	960	DO i = 1, klon
	961	! Overwrite phase partitioning in boundary layer mixed phase clouds when the
	962	! iflag_cloudth_vert=7 and specific param is activated
	963	IF (mpc_bl_points(i, k) > 0) THEN
	964	zcond(i) = MAX(0.0, qincloud_mpc(i)) * rneb(i, k)
	965	! following line is very strange and probably wrong
	966	rhcl(i, k) = (zqs(i) + zq(i)) / 2. / zqs(i)
	967	! water vapor update and partition function if thermals
	968	zq(i) = zq(i) - zcond(i)
	969	zfice(i) = zfice_th(i)
	970	ELSE
	971	! Checks on rneb, rhcl and zqn
	972	IF (rneb(i, k) <= 0.0) THEN
	973	zqn(i) = 0.0
	974	rneb(i, k) = 0.0
	975	zcond(i) = 0.0
	976	rhcl(i, k) = zq(i) / zqs(i)
	977	ELSE IF (rneb(i, k) >= 1.0) THEN
	978	zqn(i) = zq(i)
	979	rneb(i, k) = 1.0
	980	zcond(i) = MAX(0.0, zqn(i) - zqs(i))
	981	rhcl(i, k) = 1.0
	982	ELSE
	983	zcond(i) = MAX(0.0, zqn(i) - zqs(i)) * rneb(i, k)
	984	! following line is very strange and probably wrong:
	985	rhcl(i, k) = (zqs(i) + zq(i)) / 2. / zqs(i)
	986	! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
	987	IF (iflag_icefrac >= 1) THEN
	988	IF (lognormale(i)) THEN
	989	zcond(i) = zqliq(i) + zqice(i)
	990	zfice(i) = zfice_turb(i)
	991	rhcl(i, k) = zqvapcl(i) * rneb(i, k) + (zq(i) - zqn(i)) * (1. - rneb(i, k))
	992	ENDIF
[3999]	993	ENDIF
[5111]	994	ENDIF
[3999]	995
[5111]	996	! water vapor update
	997	zq(i) = zq(i) - zcond(i)
[4392]	998
[3999]	999	ENDIF
	1000
[5111]	1001	! c_iso : routine that computes in-cloud supersaturation
	1002	! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
[4226]	1003
[5111]	1004	! temperature update due to phase change
	1005	zt(i) = zt(i) + (1. - zfice(i)) * zcond(i) &
	1006	* RLVTT / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i) + zcond(i))) &
	1007	+ zfice(i) * zcond(i) * RLSTT / RCPD / (1.0 + RVTMP2 * (zq(i) + zmqc(i) + zcond(i)))
[4803]	1008	ENDDO
	1009
[5111]	1010	! If vertical heterogeneity, change volume fraction
	1011	IF (iflag_cloudth_vert >= 3) THEN
	1012	ctot_vol(1:klon, k) = min(max(ctot_vol(1:klon, k), 0.), 1.)
	1013	rneblsvol(1:klon, k) = ctot_vol(1:klon, k)
	1014	ENDIF
[4803]	1015
[5111]	1016	! ----------------------------------------------------------------
	1017	! P3> Precipitation formation
	1018	! ----------------------------------------------------------------
[4819]	1019
[5111]	1020	!================================================================
	1021	IF (ok_poprecip) THEN
[4803]	1022
[5111]	1023	DO i = 1, klon
	1024	zoliql(i) = zcond(i) * (1. - zfice(i))
	1025	zoliqi(i) = zcond(i) * zfice(i)
	1026	ENDDO
[4803]	1027
[5111]	1028	CALL poprecip_postcld(klon, dtime, paprs(:, k), paprs(:, k + 1), pplay(:, k), &
	1029	ctot_vol(:, k), ptconv(:, k), &
	1030	zt, zq, zoliql, zoliqi, zfice, &
	1031	rneb(:, k), znebprecipclr, znebprecipcld, &
	1032	zrfl, zrflclr, zrflcld, &
	1033	zifl, ziflclr, ziflcld, &
	1034	qraindiag(:, k), qsnowdiag(:, k), dqrauto(:, k), &
	1035	dqrcol(:, k), dqrmelt(:, k), dqrfreez(:, k), &
	1036	dqsauto(:, k), dqsagg(:, k), dqsrim(:, k), &
	1037	dqsmelt(:, k), dqsfreez(:, k) &
	1038	)
[3999]	1039
[5111]	1040	DO i = 1, klon
	1041	zoliq(i) = zoliql(i) + zoliqi(i)
	1042	IF (zoliq(i) > 0.) THEN
	1043	zfice(i) = zoliqi(i) / zoliq(i)
	1044	ELSE
	1045	zfice(i) = 0.0
	1046	ENDIF
[3999]	1047
[5111]	1048	! calculation of specific content of condensates seen by radiative scheme
	1049	IF (ok_radocond_snow) THEN
	1050	radocond(i, k) = zoliq(i)
	1051	radocondl(i, k) = radocond(i, k) * (1. - zfice(i))
	1052	radocondi(i, k) = radocond(i, k) * zfice(i) + qsnowdiag(i, k)
	1053	ELSE
	1054	radocond(i, k) = zoliq(i)
	1055	radocondl(i, k) = radocond(i, k) * (1. - zfice(i))
	1056	radocondi(i, k) = radocond(i, k) * zfice(i)
	1057	ENDIF
	1058	ENDDO
[3999]	1059
[5111]	1060	!================================================================
	1061	ELSE
[3999]	1062
[5111]	1063	! LTP:
	1064	IF (iflag_evap_prec >= 4) THEN
[3999]	1065
[5111]	1066	!Partitionning between precipitation coming from clouds and that coming from CS
[3999]	1067
[5111]	1068	!0) Calculate tot_zneb, total cloud fraction above the cloud
	1069	!assuming a maximum-random overlap (voir Jakob and Klein, 2000)
[3999]	1070
[5111]	1071	DO i = 1, klon
	1072	tot_znebn(i) = 1. - (1. - tot_zneb(i)) * (1 - max(rneb(i, k), zneb(i))) &
	1073	/ (1. - min(zneb(i), 1. - smallestreal))
	1074	d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
	1075	tot_zneb(i) = tot_znebn(i)
[3999]	1076
[4559]	1077
[5111]	1078	!1) Cloudy to clear air
	1079	d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i, k), znebprecipcld(i))
	1080	IF (znebprecipcld(i) > 0.) THEN
	1081	d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i) / znebprecipcld(i) * zrflcld(i)
	1082	d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i) / znebprecipcld(i) * ziflcld(i)
	1083	ELSE
	1084	d_zrfl_cld_clr(i) = 0.
	1085	d_zifl_cld_clr(i) = 0.
	1086	ENDIF
[4559]	1087
[5111]	1088	!2) Clear to cloudy air
	1089	d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i, k) - d_tot_zneb(i) - zneb(i)))
	1090	IF (znebprecipclr(i) > 0) THEN
	1091	d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i) / znebprecipclr(i) * zrflclr(i)
	1092	d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i) / znebprecipclr(i) * ziflclr(i)
	1093	ELSE
	1094	d_zrfl_clr_cld(i) = 0.
	1095	d_zifl_clr_cld(i) = 0.
	1096	ENDIF
[4559]	1097
[5111]	1098	!Update variables
	1099	znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
	1100	znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
	1101	zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
	1102	ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
	1103	zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
	1104	ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
[3999]	1105
[5111]	1106	! c_iso do the same thing for isotopes precip
	1107	ENDDO
	1108	ENDIF
[3999]	1109
	1110
[5111]	1111	! Autoconversion
	1112	! -------------------------------------------------------------------------------
[3999]	1113
	1114
[5111]	1115	! Initialisation of zoliq and radocond variables
	1116
[3999]	1117	DO i = 1, klon
[5111]	1118	zoliq(i) = zcond(i)
	1119	zoliqi(i) = zoliq(i) * zfice(i)
	1120	zoliql(i) = zoliq(i) * (1. - zfice(i))
	1121	! c_iso : initialisation of zoliq* also for isotopes
	1122	zrho(i, k) = pplay(i, k) / zt(i) / RD
	1123	zdz(i) = (paprs(i, k) - paprs(i, k + 1)) / (zrho(i, k) * RG)
	1124	iwc(i) = 0.
	1125	zneb(i) = MAX(rneb(i, k), seuil_neb)
	1126	radocond(i, k) = zoliq(i) / REAL(niter_lscp + 1)
	1127	radocondi(i, k) = zoliq(i) * zfice(i) / REAL(niter_lscp + 1)
	1128	radocondl(i, k) = zoliq(i) * (1. - zfice(i)) / REAL(niter_lscp + 1)
	1129	ENDDO
[3999]	1130
[5111]	1131	DO n = 1, niter_lscp
[3999]	1132
[5111]	1133	DO i = 1, klon
	1134	IF (rneb(i, k)>0.0) THEN
	1135	iwc(i) = zrho(i, k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
	1136	ENDIF
	1137	ENDDO
[3999]	1138
[5111]	1139	CALL fallice_velocity(klon, iwc(:), zt(:), zrho(:, k), pplay(:, k), ptconv(:, k), velo(:, k))
[4420]	1140
[5111]	1141	DO i = 1, klon
[3999]	1142
[5111]	1143	IF (rneb(i, k)>0.0) THEN
[4559]	1144
[5111]	1145	! Initialization of zrain, zsnow and zprecip:
	1146	zrain = 0.
	1147	zsnow = 0.
	1148	zprecip = 0.
	1149	! c_iso same init for isotopes. Externalisation?
[4072]	1150
[5111]	1151	IF (zneb(i)==seuil_neb) THEN
	1152	zprecip = 0.0
	1153	zsnow = 0.0
	1154	zrain = 0.0
	1155	ELSE
[3999]	1156
[5111]	1157	IF (ptconv(i, k)) THEN ! if convective point
	1158	zcl = cld_lc_con
	1159	zct = 1. / cld_tau_con
	1160	zexpo = cld_expo_con
	1161	ffallv = ffallv_con
	1162	ELSE
	1163	zcl = cld_lc_lsc
	1164	zct = 1. / cld_tau_lsc
	1165	zexpo = cld_expo_lsc
	1166	ffallv = ffallv_lsc
	1167	ENDIF
[4420]	1168
	1169
[5111]	1170	! if vertical heterogeneity is taken into account, we use
	1171	! the "true" volume fraction instead of a modified
	1172	! surface fraction (which is larger and artificially
	1173	! reduces the in-cloud water).
[4420]	1174
[5111]	1175	! Liquid water quantity to remove: zchau (Sundqvist, 1978)
	1176	! dqliq/dt=-qliq/tau(1-exp(-qcin/clw)*2)
	1177	!.........................................................
	1178	IF ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) THEN
[4420]	1179
[5111]	1180	! if vertical heterogeneity is taken into account, we use
	1181	! the "true" volume fraction instead of a modified
	1182	! surface fraction (which is larger and artificially
	1183	! reduces the in-cloud water).
	1184	effective_zneb = ctot_vol(i, k)
	1185	ELSE
	1186	effective_zneb = zneb(i)
	1187	ENDIF
[4559]	1188
[5111]	1189	IF (iflag_autoconversion == 2) THEN
	1190	! two-steps resolution with niter_lscp=1 sufficient
	1191	! we first treat the second term (with exponential) in an explicit way
	1192	! and then treat the first term (-q/tau) in an exact way
	1193	zchau = zoliql(i) * (1. - exp(-dtime / REAL(niter_lscp) * zct &
	1194	* (1. - exp(-(zoliql(i) / effective_zneb / zcl)**zexpo))))
	1195	ELSE
	1196	! old explicit resolution with subtimesteps
	1197	zchau = zct * dtime / REAL(niter_lscp) * zoliql(i) &
	1198	* (1.0 - EXP(-(zoliql(i) / effective_zneb / zcl)**zexpo))
[3999]	1199	ENDIF
[4226]	1200
[5111]	1201
	1202	! Ice water quantity to remove (Zender & Kiehl, 1997)
	1203	! dqice/dt=1/rhod(rhowice*qice)/dz
	1204	!....................................
	1205	IF (iflag_autoconversion == 2) THEN
	1206	! exact resolution, niter_lscp=1 is sufficient but works only
	1207	! with iflag_vice=0
	1208	IF (zoliqi(i) > 0.) THEN
	1209	zfroi = (zoliqi(i) - ((zoliqi(i)**(-dice_velo)) &
	1210	+ dice_velo * dtime / REAL(niter_lscp) * cice_velo / zdz(i) * ffallv)**(-1. / dice_velo))
	1211	ELSE
	1212	zfroi = 0.
	1213	ENDIF
[4559]	1214	ELSE
[5111]	1215	! old explicit resolution with subtimesteps
	1216	zfroi = dtime / REAL(niter_lscp) / zdz(i) * zoliqi(i) * velo(i, k)
[4559]	1217	ENDIF
[3999]	1218
[5111]	1219	zrain = MIN(MAX(zchau, 0.0), zoliql(i))
	1220	zsnow = MIN(MAX(zfroi, 0.0), zoliqi(i))
	1221	zprecip = MAX(zrain + zsnow, 0.0)
[4392]	1222
[5111]	1223	ENDIF
[4226]	1224
[5111]	1225	IF (iflag_autoconversion >= 1) THEN
	1226	! debugged version with phase conservation through the autoconversion process
	1227	zoliql(i) = MAX(zoliql(i) - 1. * zrain, 0.0)
	1228	zoliqi(i) = MAX(zoliqi(i) - 1. * zsnow, 0.0)
	1229	zoliq(i) = MAX(zoliq(i) - zprecip, 0.0)
	1230	ELSE
	1231	! bugged version with phase resetting
	1232	zoliql(i) = MAX(zoliq(i) * (1. - zfice(i)) - 1. * zrain, 0.0)
	1233	zoliqi(i) = MAX(zoliq(i) * zfice(i) - 1. * zsnow, 0.0)
	1234	zoliq(i) = MAX(zoliq(i) - zprecip, 0.0)
	1235	ENDIF
	1236
	1237	! c_iso: CALL isotope_conversion (for readibility) or duplicate
	1238
	1239	radocond(i, k) = radocond(i, k) + zoliq(i) / REAL(niter_lscp + 1)
	1240	radocondl(i, k) = radocondl(i, k) + zoliql(i) / REAL(niter_lscp + 1)
	1241	radocondi(i, k) = radocondi(i, k) + zoliqi(i) / REAL(niter_lscp + 1)
	1242
[4072]	1243	ENDIF ! rneb >0
[3999]	1244
[5111]	1245	ENDDO ! i = 1,klon
[3999]	1246
[5111]	1247	ENDDO ! n = 1,niter
[3999]	1248
[5111]	1249	! Precipitation flux calculation
[3999]	1250
[5111]	1251	DO i = 1, klon
[3999]	1252
[5111]	1253	IF (iflag_evap_prec>=4) THEN
	1254	ziflprev(i) = ziflcld(i)
	1255	ELSE
	1256	ziflprev(i) = zifl(i) * zneb(i)
	1257	ENDIF
[3999]	1258
[5111]	1259	IF (rneb(i, k) > 0.0) THEN
[4072]	1260
[5111]	1261	! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
	1262	! If T<0C, liquid precip are converted into ice, which leads to an increase in
	1263	! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
	1264	! taken into account through a linearization of the equations and by approximating
	1265	! the liquid precipitation process with a "threshold" process. We assume that
	1266	! condensates are not modified during this operation. Liquid precipitation is
	1267	! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
	1268	! Water vapor increases as well
	1269	! Note that compared to fisrtilp, we always assume iflag_bergeron=2
[4226]	1270
[5111]	1271	zqpreci(i) = (zcond(i) - zoliq(i)) * zfice(i)
	1272	zqprecl(i) = (zcond(i) - zoliq(i)) * (1. - zfice(i))
	1273	zcp = RCPD * (1.0 + RVTMP2 * (zq(i) + zmqc(i) + zcond(i)))
	1274	coef1 = rneb(i, k) * RLSTT / zcp * zdqsdT_raw(i)
	1275	! Computation of DT if all the liquid precip freezes
	1276	DeltaT = RLMLT * zqprecl(i) / (zcp * (1. + coef1))
	1277	! T should not exceed the freezing point
	1278	! that is Delta > RTT-zt(i)
	1279	DeltaT = max(min(RTT - zt(i), DeltaT), 0.)
	1280	zt(i) = zt(i) + DeltaT
	1281	! water vaporization due to temp. increase
	1282	Deltaq = rneb(i, k) * zdqsdT_raw(i) * DeltaT
	1283	! we add this vaporized water to the total vapor and we remove it from the precip
	1284	zq(i) = zq(i) + Deltaq
	1285	! The three "max" lines herebelow protect from rounding errors
	1286	zcond(i) = max(zcond(i) - Deltaq, 0.)
	1287	! liquid precipitation converted to ice precip
	1288	Deltaqprecl = -zcp / RLMLT * (1. + coef1) * DeltaT
	1289	zqprecl(i) = max(zqprecl(i) + Deltaqprecl, 0.)
	1290	! iced water budget
	1291	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
[4392]	1292
[5111]	1293	! c_iso : mv here condensation of isotopes + redispatchage en precip
[3999]	1294
[5111]	1295	IF (iflag_evap_prec>=4) THEN
	1296	zrflcld(i) = zrflcld(i) + zqprecl(i) &
	1297	* (paprs(i, k) - paprs(i, k + 1)) / (RG * dtime)
	1298	ziflcld(i) = ziflcld(i) + zqpreci(i) &
	1299	* (paprs(i, k) - paprs(i, k + 1)) / (RG * dtime)
	1300	znebprecipcld(i) = rneb(i, k)
	1301	zrfl(i) = zrflcld(i) + zrflclr(i)
	1302	zifl(i) = ziflcld(i) + ziflclr(i)
	1303	ELSE
	1304	zrfl(i) = zrfl(i) + zqprecl(i) &
	1305	* (paprs(i, k) - paprs(i, k + 1)) / (RG * dtime)
	1306	zifl(i) = zifl(i) + zqpreci(i) &
	1307	* (paprs(i, k) - paprs(i, k + 1)) / (RG * dtime)
	1308	ENDIF
	1309	! c_iso : same for isotopes
[3999]	1310
[5111]	1311	ENDIF ! rneb>0
[4114]	1312
[5111]	1313	ENDDO
[4114]	1314
[5111]	1315	! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
	1316	! if iflag_evap_prec>=4
	1317	IF (iflag_evap_prec>=4) THEN
	1318
	1319	DO i = 1, klon
	1320
	1321	IF ((zrflclr(i) + ziflclr(i)) > 0.) THEN
	1322	znebprecipclr(i) = min(znebprecipclr(i), max(zrflclr(i) / &
	1323	(MAX(znebprecipclr(i), seuil_neb) * rain_int_min), ziflclr(i) / (MAX(znebprecipclr(i), seuil_neb) * rain_int_min)))
[3999]	1324	ELSE
[5111]	1325	znebprecipclr(i) = 0.0
[4226]	1326	ENDIF
	1327
[5082]	1328	IF ((zrflcld(i) + ziflcld(i)) > 0.) THEN
[5111]	1329	znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i) / &
	1330	(MAX(znebprecipcld(i), seuil_neb) * rain_int_min), ziflcld(i) / (MAX(znebprecipcld(i), seuil_neb) * rain_int_min)))
[4226]	1331	ELSE
[5111]	1332	znebprecipcld(i) = 0.0
[4226]	1333	ENDIF
[5194]	1334	ENDDO
[3999]	1335
[5111]	1336	ENDIF
[3999]	1337
[5111]	1338	ENDIF ! ok_poprecip
[4910]	1339
[5111]	1340	! End of precipitation formation
	1341	! --------------------------------
[4803]	1342
[3999]	1343
[5111]	1344	! Calculation of cloud condensates amount seen by radiative scheme
	1345	!-----------------------------------------------------------------
[3999]	1346
[5111]	1347	! Calculation of the concentration of condensates seen by the radiation scheme
	1348	! for liquid phase, we take radocondl
	1349	! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
	1350	! we recalculate radocondi to account for contributions from the precipitation flux
	1351	! TODO: for the moment, we deactivate this option when ok_poprecip=.TRUE.
[4803]	1352
[5111]	1353	IF ((ok_radocond_snow) .AND. (k < klev) .AND. (.NOT. ok_poprecip)) THEN
[4114]	1354	! for the solid phase (crystals + snowflakes)
[4412]	1355	! we recalculate radocondi by summing
[4114]	1356	! the ice content calculated in the mesh
[5111]	1357	! + the contribution of the non-evaporated snowfall
[4114]	1358	! from the overlying layer
[5111]	1359	DO i = 1, klon
	1360	IF (ziflprev(i) /= 0.0) THEN
	1361	radocondi(i, k) = zoliq(i) * zfice(i) + zqpreci(i) + ziflprev(i) / zrho(i, k + 1) / velo(i, k + 1)
	1362	ELSE
	1363	radocondi(i, k) = zoliq(i) * zfice(i) + zqpreci(i)
	1364	ENDIF
	1365	radocond(i, k) = radocondl(i, k) + radocondi(i, k)
[4114]	1366	ENDDO
[5111]	1367	ENDIF
[4114]	1368
[5111]	1369	! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
	1370	DO i = 1, klon
	1371	IF (radocond(i, k) > 0.) THEN
	1372	radicefrac(i, k) = MIN(MAX(radocondi(i, k) / radocond(i, k), 0.), 1.)
[4114]	1373	ENDIF
[5111]	1374	ENDDO
[4114]	1375
[5111]	1376	! ----------------------------------------------------------------
	1377	! P4> Wet scavenging
	1378	! ----------------------------------------------------------------
[3999]	1379
[5111]	1380	!Scavenging through nucleation in the layer
[3999]	1381
[5111]	1382	DO i = 1, klon
	1383
	1384	IF(zcond(i)>zoliq(i) + 1.e-10) THEN
	1385	beta(i, k) = (zcond(i) - zoliq(i)) / zcond(i) / dtime
[3999]	1386	ELSE
[5111]	1387	beta(i, k) = 0.
[3999]	1388	ENDIF
	1389
[5111]	1390	zprec_cond(i) = MAX(zcond(i) - zoliq(i), 0.0) * (paprs(i, k) - paprs(i, k + 1)) / RG
[3999]	1391
[5111]	1392	IF (rneb(i, k)>0.0.AND.zprec_cond(i)>0.) THEN
[3999]	1393
[5111]	1394	IF (temp(i, k) >= t_glace_min) THEN
	1395	zalpha_tr = a_tr_sca(3)
	1396	ELSE
	1397	zalpha_tr = a_tr_sca(4)
	1398	ENDIF
[3999]	1399
[5111]	1400	zfrac_lessi = 1. - EXP(zalpha_tr * zprec_cond(i) / zneb(i))
	1401	frac_nucl(i, k) = 1. - zneb(i) * zfrac_lessi
[4226]	1402
[5111]	1403	! Nucleation with a factor of -1 instead of -0.5
	1404	zfrac_lessi = 1. - EXP(-zprec_cond(i) / zneb(i))
[3999]	1405
	1406	ENDIF
	1407
[5111]	1408	ENDDO
[4226]	1409
[5111]	1410	! Scavenging through impaction in the underlying layer
[3999]	1411
[5111]	1412	DO kk = k - 1, 1, -1
[3999]	1413
	1414	DO i = 1, klon
	1415
[5111]	1416	IF (rneb(i, k)>0.0.AND.zprec_cond(i)>0.) THEN
[3999]	1417
[5111]	1418	IF (temp(i, kk) >= t_glace_min) THEN
	1419	zalpha_tr = a_tr_sca(1)
	1420	ELSE
	1421	zalpha_tr = a_tr_sca(2)
	1422	ENDIF
[3999]	1423
[5111]	1424	zfrac_lessi = 1. - EXP(zalpha_tr * zprec_cond(i) / zneb(i))
	1425	frac_impa(i, kk) = 1. - zneb(i) * zfrac_lessi
[3999]	1426
[5111]	1427	ENDIF
[3999]	1428
	1429	ENDDO
	1430
[5111]	1431	ENDDO
[4226]	1432
[5111]	1433	!--save some variables for ice supersaturation
[5099]	1434
[5111]	1435	DO i = 1, klon
	1436	! for memory
	1437	rneb_seri(i, k) = rneb(i, k)
[3999]	1438
[4072]	1439	! for diagnostics
[5111]	1440	rnebclr(i, k) = 1.0 - rneb(i, k) - rnebss(i, k)
[3999]	1441
[5111]	1442	qvc(i, k) = zqs(i) * rneb(i, k)
	1443	qclr(i, k) = MAX(1.e-10, zq(i) - qvc(i, k) - qss(i, k)) !--added by OB because of pathologic cases with the lognormal
	1444	qcld(i, k) = qvc(i, k) + zcond(i)
	1445	ENDDO
	1446	!q_sat
	1447	CALL calc_qsat_ecmwf(klon, Tbef(:), qzero(:), pplay(:, k), RTT, 1, .FALSE., qsatl(:, k), zdqs(:))
	1448	CALL calc_qsat_ecmwf(klon, Tbef(:), qzero(:), pplay(:, k), RTT, 2, .FALSE., qsats(:, k), zdqs(:))
[4059]	1449
[4803]	1450
	1451
[5111]	1452	! Outputs:
	1453	!-------------------------------
	1454	! Precipitation fluxes at layer interfaces
	1455	! + precipitation fractions +
	1456	! temperature and water species tendencies
	1457	DO i = 1, klon
	1458	psfl(i, k) = zifl(i)
	1459	prfl(i, k) = zrfl(i)
	1460	pfraclr(i, k) = znebprecipclr(i)
	1461	pfracld(i, k) = znebprecipcld(i)
	1462	d_ql(i, k) = (1 - zfice(i)) * zoliq(i)
	1463	d_qi(i, k) = zfice(i) * zoliq(i)
	1464	d_q(i, k) = zq(i) - qt(i, k)
[4803]	1465	! c_iso: same for isotopes
[5111]	1466	d_t(i, k) = zt(i) - temp(i, k)
	1467	ENDDO
	1468
[4803]	1469	ENDDO
	1470
	1471
[5111]	1472	! Rain or snow at the surface (depending on the first layer temperature)
	1473	DO i = 1, klon
[3999]	1474	snow(i) = zifl(i)
	1475	rain(i) = zrfl(i)
[4392]	1476	! c_iso final output
[5111]	1477	ENDDO
[3999]	1478
[5111]	1479	IF (ncoreczq>0) THEN
	1480	WRITE(lunout, *)'WARNING : ZQ in LSCP ', ncoreczq, ' val < 1.e-15.'
	1481	ENDIF
[3999]	1482
[5111]	1483	END SUBROUTINE lscp
	1484	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4654]	1485
[4664]	1486	END MODULE lmdz_lscp

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format