Context Navigation

lmdz_lscp.f90 @ 5573

Last change on this file since 5573 was 5573, checked in by aborella, 4 months ago
Interpolation of aviation file (to be fixed) JJQ
File size: 52.6 KB

Rev	Line
[4664]	1	MODULE lmdz_lscp
[3999]	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
	7	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4380]	8	SUBROUTINE lscp(klon,klev,dtime,missing_val, &
[5396]	9	paprs, pplay, omega, temp, qt, ql_seri, qi_seri, &
	10	ptconv, ratqs, sigma_qtherm, &
[5204]	11	d_t, d_q, d_ql, d_qi, rneb, rneblsvol, &
[5007]	12	pfraclr, pfracld, &
	13	cldfraliq, sigma2_icefracturb,mean_icefracturb, &
[4530]	14	radocond, radicefrac, rain, snow, &
[3999]	15	frac_impa, frac_nucl, beta, &
[5007]	16	prfl, psfl, rhcl, qta, fraca, &
	17	tv, pspsk, tla, thl, iflag_cld_th, &
[5204]	18	iflag_ice_thermo, distcltop, temp_cltop, &
	19	tke, tke_dissip, &
	20	cell_area, &
	21	cf_seri, rvc_seri, u_seri, v_seri, &
	22	qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
[5456]	23	dcf_sub, dcf_con, dcf_mix, &
[5204]	24	dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, &
	25	dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati, &
[5573]	26	rcont_seri, flight_dist, flight_h2o, &
	27	flight_dist_read, flight_h2o_read, &
	28	aviation_lev, nleva, contfra, Tcritcont, qcritcont,&
	29	potcontfraP, potcontfraNP, dcontfra_cir, dcf_avi, &
	30	dqi_avi, dqvc_avi, cloudth_sth,cloudth_senv, &
	31	cloudth_sigmath,cloudth_sigmaenv, &
[4830]	32	qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, &
	33	dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
[4819]	34	dqsmelt, dqsfreez)
[3999]	35
	36	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	37	! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
	38	! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
	39	!--------------------------------------------------------------------------------
	40	! Date: 01/2021
	41	!--------------------------------------------------------------------------------
	42	! Aim: Large Scale Clouds and Precipitation (LSCP)
[4226]	43	!
[3999]	44	! This code is a new version of the fisrtilp.F90 routine, which is itself a
[4072]	45	! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
[3999]	46	! and 'ilp' (il pleut, L. Li)
	47	!
[4380]	48	! Compared to the original fisrtilp code, lscp
[3999]	49	! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
	50	! -> consider always precipitation thermalisation (fl_cor_ebil>0)
	51	! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
	52	! -> option "oldbug" by JYG has been removed
	53	! -> iflag_t_glace >0 always
	54	! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
	55	! -> rectangular distribution from L. Li no longer available
	56	! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
	57	!--------------------------------------------------------------------------------
	58	! References:
	59	!
[4412]	60	! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
[3999]	61	! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
	62	! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
[4412]	63	! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
[3999]	64	! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
[4412]	65	! - Touzze-Peifert Ludo, PhD thesis, p117-124
[3999]	66	! -------------------------------------------------------------------------------
	67	! Code structure:
	68	!
[4226]	69	! P0> Thermalization of the precipitation coming from the overlying layer
[3999]	70	! P1> Evaporation of the precipitation (falling from the k+1 level)
	71	! P2> Cloud formation (at the k level)
[4412]	72	! P2.A.1> With the PDFs, calculation of cloud properties using the inital
[3999]	73	! values of T and Q
	74	! P2.A.2> Coupling between condensed water and temperature
	75	! P2.A.3> Calculation of final quantities associated with cloud formation
[4412]	76	! P2.B> Release of Latent heat after cloud formation
[3999]	77	! P3> Autoconversion to precipitation (k-level)
	78	! P4> Wet scavenging
	79	!------------------------------------------------------------------------------
	80	! Some preliminary comments (JBM) :
	81	!
	82	! The cloud water that the radiation scheme sees is not the same that the cloud
	83	! water used in the physics and the dynamics
	84	!
[4412]	85	! During the autoconversion to precipitation (P3 step), radocond (cloud water used
[3999]	86	! by the radiation scheme) is calculated as an average of the water that remains
	87	! in the cloud during the precipitation and not the water remaining at the end
	88	! of the time step. The latter is used in the rest of the physics and advected
	89	! by the dynamics.
	90	!
	91	! In summary:
	92	!
	93	! Radiation:
	94	! xflwc(newmicro)+xfiwc(newmicro) =
[4412]	95	! radocond=lwcon(nc)+iwcon(nc)
[3999]	96	!
	97	! Notetheless, be aware of:
	98	!
[4412]	99	! radocond .NE. ocond(nc)
[3999]	100	! i.e.:
	101	! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
[4412]	102	! but oliq+(ocond-oliq) .EQ. ocond
[3999]	103	! (which is not trivial)
	104	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4226]	105	!
[4654]	106
	107	! USE de modules contenant des fonctions.
[4651]	108	USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
[5007]	109	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
[5413]	110	USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
[5241]	111	USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
[5413]	112	USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
	113	USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
[3999]	114
[4664]	115	! Use du module lmdz_lscp_ini contenant les constantes
[5204]	116	USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
[5413]	117	USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
	118	USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
	119	USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
	120	USE lmdz_lscp_ini, ONLY : t_glace_min, min_frac_th_cld
	121	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
[5412]	122	USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
[5211]	123	USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
[5452]	124	USE lmdz_lscp_ini, ONLY : ok_plane_contrail
[5204]	125
[5573]	126	! aviation module
	127	USE mod_phys_lmdz_para, ONLY : is_omp_master
	128	USE lmdz_aviation, ONLY : vertical_interpolation_aviation
	129
[3999]	130	IMPLICIT NONE
	131
	132	!===============================================================================
[4380]	133	! VARIABLES DECLARATION
[3999]	134	!===============================================================================
	135
	136	! INPUT VARIABLES:
	137	!-----------------
	138
[4380]	139	INTEGER, INTENT(IN) :: klon,klev ! number of horizontal grid points and vertical levels
[3999]	140	REAL, INTENT(IN) :: dtime ! time step [s]
[4380]	141	REAL, INTENT(IN) :: missing_val ! missing value for output
[4059]	142
[3999]	143	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa]
	144	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa]
[4654]	145	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K)
[5383]	146	REAL, DIMENSION(klon,klev), INTENT(IN) :: omega ! vertical velocity [Pa/s]
[4686]	147	REAL, DIMENSION(klon,klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg]
[5396]	148	REAL, DIMENSION(klon,klev), INTENT(IN) :: ql_seri ! liquid specific content [kg/kg]
	149	REAL, DIMENSION(klon,klev), INTENT(IN) :: qi_seri ! ice specific content [kg/kg]
[3999]	150	INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds
	151	INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
[4226]	152	! CR: if iflag_ice_thermo=2, only convection is active
[3999]	153	LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
	154
	155	!Inputs associated with thermal plumes
[4226]	156
[5007]	157	REAL, DIMENSION(klon,klev), INTENT(IN) :: tv ! virtual potential temperature [K]
	158	REAL, DIMENSION(klon,klev), INTENT(IN) :: qta ! specific humidity within thermals [kg/kg]
	159	REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-]
	160	REAL, DIMENSION(klon,klev), INTENT(IN) :: pspsk ! exner potential (p/100000)**(R/cp)
	161	REAL, DIMENSION(klon,klev), INTENT(IN) :: tla ! liquid temperature within thermals [K]
[5413]	162	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke !--turbulent kinetic energy [m2/s2]
	163	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip !--TKE dissipation [m2/s3]
[3999]	164
	165	! INPUT/OUTPUT variables
	166	!------------------------
[4562]	167
[5208]	168	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: thl ! liquid potential temperature [K]
	169	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs,sigma_qtherm ! function of pressure that sets the large-scale
[4059]	170
[5208]	171
[5204]	172	! INPUT/OUTPUT condensation and ice supersaturation
	173	!--------------------------------------------------
	174	REAL, DIMENSION(klon,klev), INTENT(INOUT):: cf_seri ! cloud fraction [-]
	175	REAL, DIMENSION(klon,klev), INTENT(INOUT):: rvc_seri ! cloudy water vapor to total water vapor ratio [-]
	176	REAL, DIMENSION(klon,klev), INTENT(IN) :: u_seri ! eastward wind [m/s]
	177	REAL, DIMENSION(klon,klev), INTENT(IN) :: v_seri ! northward wind [m/s]
	178	REAL, DIMENSION(klon), INTENT(IN) :: cell_area ! area of each cell [m2]
[4380]	179
[5204]	180	! INPUT/OUTPUT aviation
	181	!--------------------------------------------------
[5452]	182	REAL, DIMENSION(klon,klev), INTENT(INOUT):: rcont_seri ! ratio of contrails fraction to total cloud fraction [-]
[5573]	183	REAL, DIMENSION(klon,klev,1), INTENT(OUT) :: flight_dist ! aviation distance flown within the mesh [m/s/mesh]
	184	REAL, DIMENSION(klon,klev,1), INTENT(OUT) :: flight_h2o ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
	185	! Read from file - nleva vertical levels - define in state_var
	186	INTEGER, INTENT(IN) :: nleva
	187	REAL, DIMENSION(nleva), INTENT(IN) :: aviation_lev ! vertical levels [km]
	188	REAL, DIMENSION(klon,nleva,1), INTENT(IN) :: flight_dist_read ! aviation distance flown within the mesh [m/s/mesh]
	189	REAL, DIMENSION(klon,nleva,1), INTENT(IN) :: flight_h2o_read ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
[3999]	190
	191	! OUTPUT variables
	192	!-----------------
	193
[5204]	194	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K]
	195	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg]
	196	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg]
	197	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg]
	198	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-]
	199	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-]
	200	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-]
	201	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-]
[5007]	202	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cldfraliq ! liquid fraction of cloud [-]
	203	REAL, DIMENSION(klon,klev), INTENT(OUT) :: sigma2_icefracturb ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
	204	REAL, DIMENSION(klon,klev), INTENT(OUT) :: mean_icefracturb ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
[5204]	205	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg]
	206	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme
	207	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-]
	208	REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2]
	209	REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2]
	210	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2]
	211	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2]
	212	REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m]
	213	REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K]
	214	REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water
[4686]	215
[5007]	216	! fraction of aerosol scavenging through impaction and nucleation (for on-line)
[3999]	217
[5007]	218	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-]
	219	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-]
[3999]	220
[5204]	221	! for condensation and ice supersaturation
[3999]	222
[5204]	223	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsub !--specific total water content in sub-saturated clear sky region [kg/kg]
	224	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qissr !--specific total water content in supersat region [kg/kg]
	225	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld !--specific total water content in cloudy region [kg/kg]
	226	REAL, DIMENSION(klon,klev), INTENT(OUT) :: subfra !--mesh fraction of subsaturated clear sky [-]
	227	REAL, DIMENSION(klon,klev), INTENT(OUT) :: issrfra !--mesh fraction of ISSR [-]
	228	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_cond !--coefficient governing the ice nucleation RHi threshold [-]
	229	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_sub !--cloud fraction tendency because of sublimation [s-1]
	230	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_con !--cloud fraction tendency because of condensation [s-1]
	231	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_mix !--cloud fraction tendency because of cloud mixing [s-1]
	232	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_adj !--specific ice content tendency because of temperature adjustment [kg/kg/s]
	233	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_sub !--specific ice content tendency because of sublimation [kg/kg/s]
	234	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_con !--specific ice content tendency because of condensation [kg/kg/s]
	235	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_mix !--specific ice content tendency because of cloud mixing [kg/kg/s]
	236	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_adj !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
	237	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_sub !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
	238	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_con !--specific cloud water vapor tendency because of condensation [kg/kg/s]
	239	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_mix !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
	240	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl !--saturation specific humidity wrt liquid [kg/kg]
	241	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsati !--saturation specific humidity wrt ice [kg/kg]
	242
	243	! for contrails and aviation
	244
[5456]	245	REAL, DIMENSION(klon,klev), INTENT(OUT) :: contfra !--linear contrail fraction [-]
	246	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcritcont !--critical temperature for contrail formation [K]
	247	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcritcont !--critical specific humidity for contrail formation [kg/kg]
	248	REAL, DIMENSION(klon,klev), INTENT(OUT) :: potcontfraP !--potential persistent contrail fraction [-]
	249	REAL, DIMENSION(klon,klev), INTENT(OUT) :: potcontfraNP !--potential non-persistent contrail fraction [-]
	250	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcontfra_cir !--linear contrail fraction to cirrus cloud fraction tendency [s-1]
[5204]	251	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dcf_avi !--cloud fraction tendency because of aviation [s-1]
	252	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqi_avi !--specific ice content tendency because of aviation [kg/kg/s]
	253	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqvc_avi !--specific cloud water vapor tendency because of aviation [kg/kg/s]
	254
	255
[4803]	256	! for POPRECIP
[4686]	257
[4830]	258	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg]
	259	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg]
[4819]	260	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqreva !--rain tendendy due to evaporation [kg/kg/s]
	261	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
	262	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
	263	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
	264	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
	265	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
	266	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsrim !--snow tendency due to riming [kg/kg/s]
	267	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsmelt !--snow tendency due to melting [kg/kg/s]
	268	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrmelt !--rain tendency due to melting [kg/kg/s]
	269	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsfreez !--snow tendency due to freezing [kg/kg/s]
	270	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrfreez !--rain tendency due to freezing [kg/kg/s]
[4686]	271
[5204]	272	! for thermals
[4803]	273
[5204]	274	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth !--mean saturation deficit in thermals
	275	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv !--mean saturation deficit in environment
	276	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath !--std of saturation deficit in thermals
	277	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv !--std of saturation deficit in environment
	278
	279
[3999]	280	! LOCAL VARIABLES:
	281	!----------------
[5383]	282	REAL,DIMENSION(klon) :: qice_ini
[4654]	283	REAL, DIMENSION(klon,klev) :: ctot
	284	REAL, DIMENSION(klon,klev) :: ctot_vol
[5450]	285	REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi
[5413]	286	REAL :: zdelta
[4654]	287	REAL, DIMENSION(klon) :: zdqsdT_raw
[5007]	288	REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T
	289	REAL, DIMENSION(klon) :: Tbef,qlbef,DT ! temperature, humidity and temp. variation during lognormal iteration
[4654]	290	REAL :: num,denom
	291	REAL :: cste
[4910]	292	REAL, DIMENSION(klon) :: zfice_th
[4654]	293	REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
[5413]	294	REAL, DIMENSION(klon) :: zrfl, zifl
[5450]	295	REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn, zqnl
[4654]	296	REAL, DIMENSION(klon) :: zoliql, zoliqi
	297	REAL, DIMENSION(klon) :: zt
[5450]	298	REAL, DIMENSION(klon) :: zfice,zneb, znebl
[4654]	299	REAL, DIMENSION(klon) :: dzfice
[5007]	300	REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
[4654]	301	REAL, DIMENSION(klon) :: qtot, qzero
[3999]	302	! Variables precipitation energy conservation
[4654]	303	REAL, DIMENSION(klon) :: zmqc
	304	REAL :: zalpha_tr
	305	REAL :: zfrac_lessi
	306	REAL, DIMENSION(klon) :: zprec_cond
	307	REAL, DIMENSION(klon) :: zlh_solid
[4803]	308	REAL, DIMENSION(klon) :: ztupnew
[4913]	309	REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
[4654]	310	REAL, DIMENSION(klon) :: zrflclr, zrflcld
	311	REAL, DIMENSION(klon) :: ziflclr, ziflcld
[5411]	312	REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
[5413]	313	REAL, DIMENSION(klon) :: tot_zneb
[4654]	314	REAL :: qlmpc, qimpc, rnebmpc
[5007]	315	REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
	316	REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl ! for icefrac_lscp_turb
[5413]	317
	318	! for quantity of condensates seen by radiation
	319	REAL, DIMENSION(klon) :: zradocond, zradoice
	320	REAL, DIMENSION(klon) :: zrho_up, zvelo_up
[5204]	321
	322	! for condensation and ice supersaturation
[5450]	323	REAL, DIMENSION(klon) :: qvc, qvcl, shear
[5204]	324	REAL :: delta_z
[5452]	325	! for contrails
	326	!--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrail)
[5204]	327	! Constants used for calculating ratios that are advected (using a parent-child
	328	! formalism). This is not done in the dynamical core because at this moment,
	329	! only isotopes can use this parent-child formalism. Note that the two constants
	330	! are the same as the one use in the dynamical core, being also defined in
	331	! dyn3d_common/infotrac.F90
	332	REAL :: min_qParent, min_ratio
[3999]	333
[5413]	334	INTEGER i, k, kk, iter
[4654]	335	INTEGER, DIMENSION(klon) :: n_i
[4226]	336	INTEGER ncoreczq
[4654]	337	INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
[4803]	338	LOGICAL iftop
[3999]	339
[4654]	340	LOGICAL, DIMENSION(klon) :: lognormale
	341	LOGICAL, DIMENSION(klon) :: keepgoing
[3999]	342
	343	CHARACTER (len = 20) :: modname = 'lscp'
	344	CHARACTER (len = 80) :: abort_message
	345
	346
	347	!===============================================================================
	348	! INITIALISATION
	349	!===============================================================================
	350
[4226]	351	! Few initial checks
[3999]	352
[4654]	353
[3999]	354	IF (iflag_fisrtilp_qsat .LT. 0) THEN
[4226]	355	abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
	356	CALL abort_physic(modname,abort_message,1)
[3999]	357	ENDIF
	358
	359	! Few initialisations
	360
	361	ctot_vol(1:klon,1:klev)=0.0
	362	rneblsvol(1:klon,1:klev)=0.0
[5411]	363	znebprecip(:)=0.0
[4226]	364	znebprecipclr(:)=0.0
	365	znebprecipcld(:)=0.0
[3999]	366	mpc_bl_points(:,:)=0
	367
	368	IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
	369
	370	! AA for 'safety' reasons
	371	zalpha_tr = 0.
	372	zfrac_lessi = 0.
[4380]	373	beta(:,:)= 0.
[3999]	374
[4380]	375	! Initialisation of variables:
	376
[3999]	377	prfl(:,:) = 0.0
	378	psfl(:,:) = 0.0
	379	d_t(:,:) = 0.0
	380	d_q(:,:) = 0.0
	381	d_ql(:,:) = 0.0
	382	d_qi(:,:) = 0.0
	383	rneb(:,:) = 0.0
[4530]	384	pfraclr(:,:)=0.0
	385	pfracld(:,:)=0.0
[5007]	386	cldfraliq(:,:)=0.
	387	sigma2_icefracturb(:,:)=0.
	388	mean_icefracturb(:,:)=0.
[4412]	389	radocond(:,:) = 0.0
[3999]	390	radicefrac(:,:) = 0.0
[4226]	391	frac_nucl(:,:) = 1.0
	392	frac_impa(:,:) = 1.0
[3999]	393	rain(:) = 0.0
	394	snow(:) = 0.0
[5450]	395	zfice(:)=1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
[4114]	396	dzfice(:)=0.0
[5007]	397	zfice_turb(:)=0.0
	398	dzfice_turb(:)=0.0
[3999]	399	zrfl(:) = 0.0
	400	zifl(:) = 0.0
	401	zneb(:) = seuil_neb
[4226]	402	zrflclr(:) = 0.0
	403	ziflclr(:) = 0.0
	404	zrflcld(:) = 0.0
	405	ziflcld(:) = 0.0
	406	tot_zneb(:) = 0.0
	407	qzero(:) = 0.0
[5007]	408	zdistcltop(:)=0.0
	409	ztemp_cltop(:) = 0.0
[4803]	410	ztupnew(:)=0.0
[5204]	411
[4654]	412	distcltop(:,:)=0.
	413	temp_cltop(:,:)=0.
[5007]	414
[5204]	415	!--Ice supersaturation
	416	gamma_cond(:,:) = 1.
	417	qissr(:,:) = 0.
	418	issrfra(:,:) = 0.
	419	dcf_sub(:,:) = 0.
	420	dcf_con(:,:) = 0.
	421	dcf_mix(:,:) = 0.
	422	dqi_adj(:,:) = 0.
	423	dqi_sub(:,:) = 0.
	424	dqi_con(:,:) = 0.
	425	dqi_mix(:,:) = 0.
	426	dqvc_adj(:,:) = 0.
	427	dqvc_sub(:,:) = 0.
	428	dqvc_con(:,:) = 0.
	429	dqvc_mix(:,:) = 0.
[5456]	430	contfra(:,:) = 0.
[5452]	431	Tcritcont(:,:) = missing_val
	432	qcritcont(:,:) = missing_val
[5455]	433	potcontfraP(:,:)= 0.
	434	potcontfraNP(:,:)= 0.
[5456]	435	dcontfra_cir(:,:)= 0.
[5204]	436	dcf_avi(:,:) = 0.
	437	dqi_avi(:,:) = 0.
	438	dqvc_avi(:,:) = 0.
	439	qvc(:) = 0.
	440	shear(:) = 0.
	441	min_qParent = 1.e-30
	442	min_ratio = 1.e-16
	443
[4803]	444	!-- poprecip
[4830]	445	qraindiag(:,:)= 0.
	446	qsnowdiag(:,:)= 0.
[4819]	447	dqreva(:,:) = 0.
	448	dqrauto(:,:) = 0.
	449	dqrmelt(:,:) = 0.
	450	dqrfreez(:,:) = 0.
	451	dqrcol(:,:) = 0.
	452	dqssub(:,:) = 0.
	453	dqsauto(:,:) = 0.
	454	dqsrim(:,:) = 0.
	455	dqsagg(:,:) = 0.
	456	dqsfreez(:,:) = 0.
	457	dqsmelt(:,:) = 0.
[4913]	458	zqupnew(:) = 0.
	459	zqvapclr(:) = 0.
[3999]	460
[4803]	461
	462
[4392]	463	!c_iso: variable initialisation for iso
	464
[5573]	465	IF ( ok_plane_contrail ) THEN
	466	! Vertical interpolation is done at each physical timestep
	467	!
	468	IF (is_omp_master) CALL vertical_interpolation_aviation(klon, klev, paprs, pplay, flight_dist_read, &
	469	flight_h2o_read, aviation_lev, nleva, flight_dist, flight_h2o)
	470	ENDIF
[4392]	471
[3999]	472	!===============================================================================
	473	! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
	474	!===============================================================================
	475
	476	ncoreczq=0
	477
	478	DO k = klev, 1, -1
	479
[5408]	480	qice_ini = qi_seri(:,k)
	481
[4803]	482	IF (k.LE.klev-1) THEN
	483	iftop=.false.
	484	ELSE
	485	iftop=.true.
	486	ENDIF
	487
[3999]	488	! Initialisation temperature and specific humidity
[5007]	489	! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
	490	! at the end of the klon loop, a temperature incremtent d_t due to all processes
	491	! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
	492	! d_t = temperature tendency due to lscp
	493	! The temperature of the overlying layer is updated here because needed for thermalization
[3999]	494	DO i = 1, klon
[4654]	495	zt(i)=temp(i,k)
[4686]	496	zq(i)=qt(i,k)
[4803]	497	IF (.not. iftop) THEN
[4913]	498	ztupnew(i) = temp(i,k+1) + d_t(i,k+1)
	499	zqupnew(i) = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
	500	!--zqs(i) is the saturation specific humidity in the layer above
	501	zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
[4803]	502	ENDIF
[4392]	503	!c_iso init of iso
[3999]	504	ENDDO
[5408]	505
	506	! --------------------------------------------------------------------
	507	! P1> Precipitation processes, before cloud formation:
	508	! Thermalization of precipitation falling from the overlying layer AND
	509	! Precipitation evaporation/sublimation/melting
	510	!---------------------------------------------------------------------
[4803]	511
	512	!================================================================
	513	! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
[5431]	514	IF ( ok_poprecip ) THEN
[4803]	515
[5408]	516	CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	517	zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
	518	zqvapclr, zqupnew, &
[5431]	519	cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
[5408]	520	zrfl, zrflclr, zrflcld, &
	521	zifl, ziflclr, ziflcld, &
	522	dqreva(:,k), dqssub(:,k) &
	523	)
[4803]	524
	525	!================================================================
	526	ELSE
	527
[5408]	528	CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
[5411]	529	zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
[5408]	530	zrfl, zrflclr, zrflcld, &
	531	zifl, ziflclr, ziflcld, &
	532	dqreva(:,k), dqssub(:,k) &
	533	)
[3999]	534
[4803]	535	ENDIF ! (ok_poprecip)
[3999]	536
	537	! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
	538	!-------------------------------------------------------
[4226]	539
[4072]	540	qtot(:)=zq(:)+zmqc(:)
[5383]	541	CALL calc_qsat_ecmwf(klon,zt,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
[4072]	542	DO i = 1, klon
[3999]	543	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
[4059]	544	zdqsdT_raw(i) = zdqs(i)RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTT*zdelta)
[3999]	545	IF (zq(i) .LT. 1.e-15) THEN
	546	ncoreczq=ncoreczq+1
	547	zq(i)=1.e-15
	548	ENDIF
[4392]	549	! c_iso: do something similar for isotopes
[3999]	550
[4072]	551	ENDDO
[3999]	552
	553	! --------------------------------------------------------------------
[5408]	554	! P2> Cloud formation
[3999]	555	!---------------------------------------------------------------------
	556	!
	557	! Unlike fisrtilp, we always assume a 'fractional cloud' approach
	558	! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
	559	! is prescribed and depends on large scale variables and boundary layer
	560	! properties)
	561	! The decrease in condensed part due tu latent heating is taken into
	562	! account
	563	! -------------------------------------------------------------------
	564
[5408]	565	! P2.1> With the PDFs (log-normal, bigaussian)
[4424]	566	! cloud properties calculation with the initial values of t and q
[3999]	567	! ----------------------------------------------------------------
	568
	569	! initialise gammasat and qincloud_mpc
	570	gammasat(:)=1.
	571	qincloud_mpc(:)=0.
	572
	573	IF (iflag_cld_th.GE.5) THEN
[4910]	574	! Cloud cover and content in meshes affected by shallow convection,
	575	! are retrieved from a bi-gaussian distribution of the saturation deficit
	576	! following Jam et al. 2013
[3999]	577
[4910]	578	IF (iflag_cloudth_vert.LE.2) THEN
	579	! Old version of Arnaud Jam
[3999]	580
[4910]	581	CALL cloudth(klon,klev,k,tv, &
	582	zq,qta,fraca, &
	583	qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
[4654]	584	ratqs,zqs,temp, &
[4651]	585	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	586
[4651]	587
[3999]	588	ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
[4910]	589	! Default version of Arnaud Jam
	590
	591	CALL cloudth_v3(klon,klev,k,tv, &
	592	zq,qta,fraca, &
	593	qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
[5208]	594	ratqs,sigma_qtherm,zqs,temp, &
[4651]	595	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	596
	597
	598	ELSEIF (iflag_cloudth_vert.EQ.6) THEN
[4910]	599	! Jean Jouhaud's version, with specific separation between surface and volume
	600	! cloud fraction Decembre 2018
[3999]	601
[4910]	602	CALL cloudth_v6(klon,klev,k,tv, &
	603	zq,qta,fraca, &
	604	qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
[4654]	605	ratqs,zqs,temp, &
[4651]	606	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	607
[4910]	608	ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
	609	! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
[5007]	610	! for boundary-layer mixed phase clouds
[4910]	611	CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
	612	pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
	613	ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
	614	cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
	615
[3999]	616	ENDIF
	617
	618
	619	DO i=1,klon
	620	rneb(i,k)=ctot(i,k)
	621	rneblsvol(i,k)=ctot_vol(i,k)
	622	zqn(i)=qcloud(i)
[5204]	623	!--AB grid-mean vapor in the cloud - we assume saturation adjustment
	624	qvc(i) = rneb(i,k) * zqs(i)
[3999]	625	ENDDO
	626
	627	ENDIF
	628
	629	IF (iflag_cld_th .LE. 4) THEN
	630
	631	! lognormal
[5007]	632	lognormale(:) = .TRUE.
[3999]	633
	634	ELSEIF (iflag_cld_th .GE. 6) THEN
	635
	636	! lognormal distribution when no thermals
[5007]	637	lognormale(:) = fraca(:,k) < min_frac_th_cld
[3999]	638
	639	ELSE
[4226]	640	! When iflag_cld_th=5, we always assume
[3999]	641	! bi-gaussian distribution
[5007]	642	lognormale(:) = .FALSE.
[3999]	643
	644	ENDIF
	645
	646	DT(:) = 0.
	647	n_i(:)=0
	648	Tbef(:)=zt(:)
	649	qlbef(:)=0.
	650
	651	! Treatment of non-boundary layer clouds (lognormale)
[5383]	652	! We iterate here to take into account the change in qsat(T) and ice fraction
	653	! during the condensation process
	654	! the increment in temperature is calculated using the first principle of
	655	! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
	656	! and a clear sky fraction)
	657	! note that we assume that the vapor in the cloud is at saturation for this calculation
[4226]	658
[3999]	659	DO iter=1,iflag_fisrtilp_qsat+1
[4226]	660
[5204]	661	keepgoing(:) = .FALSE.
	662
[3999]	663	DO i=1,klon
	664
[4424]	665	! keepgoing = .true. while convergence is not satisfied
[4226]	666
[5204]	667	IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
[4226]	668
[3999]	669	! if not convergence:
[5204]	670	! we calculate a new iteration
	671	keepgoing(i) = .TRUE.
[3999]	672
	673	! P2.2.1> cloud fraction and condensed water mass calculation
	674	! Calculated variables:
	675	! rneb : cloud fraction
	676	! zqn : total water within the cloud
	677	! zcond: mean condensed water within the mesh
	678	! rhcl: clear-sky relative humidity
	679	!---------------------------------------------------------------
	680
[4424]	681	! new temperature that only serves in the iteration process:
[3999]	682	Tbef(i)=Tbef(i)+DT(i)
	683
	684	! Rneb, qzn and zcond for lognormal PDFs
[4072]	685	qtot(i)=zq(i)+zmqc(i)
[4226]	686
[4072]	687	ENDIF
[3999]	688
[4072]	689	ENDDO
	690
[4380]	691	! Calculation of saturation specific humidity and ice fraction
[5383]	692	CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
[5450]	693
	694	IF (iflag_icefrac .GE. 3) THEN
	695	! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
	696	! and cloud condensed water content. idea following Dietlitcher et al. 2018, GMD
	697	! we make this option works only for the physically-based and tke-depdenent param (iflag_icefrac>=1)
	698	DO i=1,klon
	699	CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,1,.false.,zqsl,zdqsl)
	700	CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,2,.false.,zqsi,zdqsi)
	701	zqs(i)=zfice(i)zqsi(i)+(1.-zfice(i))zqsl(i)
	702	zdqs(i)=zfice(i)zdqsi(i)+zqsi(i)dzfice(i)+(1.-zfice(i))zdqsl(i)-zqsl(i)dzfice(i)
	703	ENDDO
	704	ENDIF
	705
[5383]	706	CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
[4072]	707	! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
	708	zdqs(:) = gammasat(:)zdqs(:)+zqs(:)dgammasatdt(:)
	709
[5450]	710	! Cloud condensation based on subgrid pdf
[5204]	711	!--AB Activates a condensation scheme that allows for
	712	!--ice supersaturation and contrails evolution from aviation
	713	IF (ok_ice_supersat) THEN
	714
	715	!--Calculate the shear value (input for condensation and ice supersat)
	716	DO i = 1, klon
	717	!--Cell thickness [m]
	718	delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
	719	IF ( iftop ) THEN
	720	! top
	721	shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
	722	+ ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
	723	ELSEIF ( k .EQ. 1 ) THEN
	724	! surface
	725	shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
	726	+ ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
	727	ELSE
	728	! other layers
	729	shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
	730	- (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
	731	+ ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
	732	- (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
	733	ENDIF
	734	ENDDO
	735
	736	!---------------------------------------------
	737	!-- CONDENSATION AND ICE SUPERSATURATION --
	738	!---------------------------------------------
	739
	740	CALL condensation_ice_supersat( &
	741	klon, dtime, missing_val, &
	742	pplay(:,k), paprs(:,k), paprs(:,k+1), &
[5396]	743	cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
[5406]	744	shear, tke_dissip(:,k), cell_area, &
	745	Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
	746	rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
[5204]	747	dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
	748	dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
	749	dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
[5573]	750	rcont_seri(:,k), flight_dist(:,k,1), flight_h2o(:,k,1), contfra(:,k), &
[5452]	751	Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
[5456]	752	dcontfra_cir(:,k), dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
[5204]	753
	754
[5211]	755	ELSE
	756	!--generalised lognormal condensation scheme (Bony and Emanuel 2001)
	757
[5241]	758	CALL condensation_lognormal( &
	759	klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
	760	keepgoing, rneb(:,k), zqn, qvc)
[5204]	761
	762
	763	ENDIF ! .NOT. ok_ice_supersat
[4226]	764
[5450]	765	! cloud phase determination
	766	IF (iflag_icefrac .GE. 2) THEN
	767	! phase partitioning depending on temperature. activates here in the iteration process if iflag_icefrac > 2
	768	CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
	769	rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
	770	ELSE
	771	! phase partitioning depending on temperature and eventually distance to cloud top
	772	IF (iflag_t_glace.GE.4) THEN
	773	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
	774	CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
	775	ENDIF
	776	CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
	777	ENDIF
	778
	779
[5204]	780	DO i=1,klon
	781	IF (keepgoing(i)) THEN
[3999]	782
	783	! If vertical heterogeneity, change fraction by volume as well
	784	IF (iflag_cloudth_vert.GE.3) THEN
	785	ctot_vol(i,k)=rneb(i,k)
	786	rneblsvol(i,k)=ctot_vol(i,k)
	787	ENDIF
	788
	789
[4072]	790	! P2.2.2> Approximative calculation of temperature variation DT
	791	! due to condensation.
	792	! Calculated variables:
	793	! dT : temperature change due to condensation
	794	!---------------------------------------------------------------
[3999]	795
	796
	797	IF (zfice(i).LT.1) THEN
	798	cste=RLVTT
	799	ELSE
	800	cste=RLSTT
	801	ENDIF
[5007]	802
[5204]	803	IF ( ok_unadjusted_clouds ) THEN
	804	!--AB We relax the saturation adjustment assumption
	805	!-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
	806	IF ( rneb(i,k) .GT. eps ) THEN
	807	qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
	808	ELSE
	809	qlbef(i) = 0.
	810	ENDIF
	811	ELSE
	812	qlbef(i)=max(0.,zqn(i)-zqs(i))
	813	ENDIF
	814
[3999]	815	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
[4226]	816	+zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
[3999]	817	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
	818	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
[4226]	819	qlbef(i)dzfice(i)
[4424]	820	! here we update a provisory temperature variable that only serves in the iteration
	821	! process
[3999]	822	DT(i)=num/denom
	823	n_i(i)=n_i(i)+1
	824
[4424]	825	ENDIF ! end keepgoing
[3999]	826
	827	ENDDO ! end loop on i
	828
	829	ENDDO ! iter=1,iflag_fisrtilp_qsat+1
	830
[5408]	831	! P2.2> Final quantities calculation
[3999]	832	! Calculated variables:
	833	! rneb : cloud fraction
	834	! zcond: mean condensed water in the mesh
	835	! zqn : mean water vapor in the mesh
[4562]	836	! zfice: ice fraction in clouds
[3999]	837	! zt : temperature
	838	! rhcl : clear-sky relative humidity
	839	! ----------------------------------------------------------------
	840
	841
[5450]	842	! Cloud phase final determination
	843	!--------------------------------
[4562]	844	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
	845	IF (iflag_t_glace.GE.4) THEN
[5007]	846	CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
	847	distcltop(:,k)=zdistcltop(:)
	848	temp_cltop(:,k)=ztemp_cltop(:)
	849	ENDIF
[5450]	850	! Partition function depending on temperature for all clouds (shallow convective and stratiform)
[5204]	851	CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
[4562]	852
[5450]	853	! Partition function depending on tke for non shallow-convective clouds, erase previous estimation
[5007]	854	IF (iflag_icefrac .GE. 1) THEN
[5383]	855	CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
[5007]	856	rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
	857	ENDIF
	858
[5450]	859	! Water vapor update, subsequent latent heat exchange for each cloud type
	860	!------------------------------------------------------------------------
[3999]	861	DO i=1, klon
[5007]	862	! Overwrite phase partitioning in boundary layer mixed phase clouds when the
	863	! iflag_cloudth_vert=7 and specific param is activated
[3999]	864	IF (mpc_bl_points(i,k) .GT. 0) THEN
	865	zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
	866	! following line is very strange and probably wrong
	867	rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
	868	! water vapor update and partition function if thermals
	869	zq(i) = zq(i) - zcond(i)
[4910]	870	zfice(i)=zfice_th(i)
[3999]	871	ELSE
	872	! Checks on rneb, rhcl and zqn
	873	IF (rneb(i,k) .LE. 0.0) THEN
	874	zqn(i) = 0.0
	875	rneb(i,k) = 0.0
	876	zcond(i) = 0.0
	877	rhcl(i,k)=zq(i)/zqs(i)
	878	ELSE IF (rneb(i,k) .GE. 1.0) THEN
	879	zqn(i) = zq(i)
	880	rneb(i,k) = 1.0
[5204]	881	IF ( ok_unadjusted_clouds ) THEN
	882	!--AB We relax the saturation adjustment assumption
	883	!-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
	884	zcond(i) = MAX(0., zqn(i) - qvc(i))
	885	ELSE
	886	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
	887	ENDIF
[3999]	888	rhcl(i,k)=1.0
	889	ELSE
[5204]	890	IF ( ok_unadjusted_clouds ) THEN
	891	!--AB We relax the saturation adjustment assumption
	892	!-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
	893	zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
	894	ELSE
	895	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
	896	ENDIF
[3999]	897	! following line is very strange and probably wrong:
	898	rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
[5007]	899	! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
	900	IF (iflag_icefrac .GE. 1) THEN
	901	IF (lognormale(i)) THEN
	902	zcond(i) = zqliq(i) + zqice(i)
[5450]	903	zfice(i) = zfice_turb(i)
[5007]	904	rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
	905	ENDIF
	906	ENDIF
[3999]	907	ENDIF
	908
[4226]	909	! water vapor update
	910	zq(i) = zq(i) - zcond(i)
[3999]	911
	912	ENDIF
[5450]	913
	914
[3999]	915	! temperature update due to phase change
	916	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
	917	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
	918	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
	919	ENDDO
	920
	921	! If vertical heterogeneity, change volume fraction
	922	IF (iflag_cloudth_vert .GE. 3) THEN
	923	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
	924	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
	925	ENDIF
	926
[5204]	927
[3999]	928	! ----------------------------------------------------------------
[5408]	929	! P3> Precipitation processes, after cloud formation
	930	! - precipitation formation, melting/freezing
[3999]	931	! ----------------------------------------------------------------
[4226]	932
[5413]	933	! Initiate the quantity of liquid and solid condensates
	934	! Note that in the following, zcond is the total amount of condensates
	935	! including newly formed precipitations (i.e., condensates formed by the
	936	! condensation process), while zoliq is the total amount of condensates in
	937	! the cloud (i.e., on which precipitation processes have applied)
	938	DO i = 1, klon
	939	zoliq(i) = zcond(i)
	940	zoliql(i) = zcond(i) * ( 1. - zfice(i) )
	941	zoliqi(i) = zcond(i) * zfice(i)
	942	ENDDO
	943
[4803]	944	!================================================================
[5408]	945	! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
[4803]	946	IF (ok_poprecip) THEN
	947
	948	CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	949	ctot_vol(:,k), ptconv(:,k), &
	950	zt, zq, zoliql, zoliqi, zfice, &
	951	rneb(:,k), znebprecipclr, znebprecipcld, &
	952	zrfl, zrflclr, zrflcld, &
	953	zifl, ziflclr, ziflcld, &
[4830]	954	qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
[4819]	955	dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
	956	dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
	957	dqsmelt(:,k), dqsfreez(:,k) &
[4803]	958	)
[5431]	959	DO i = 1, klon
	960	zoliq(i) = zoliql(i) + zoliqi(i)
	961	ENDDO
[4803]	962
	963	!================================================================
	964	ELSE
	965
[5413]	966	CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	967	ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
	968	zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
	969	rneb(:,k), znebprecipclr, znebprecipcld, &
	970	zneb, tot_zneb, zrho_up, zvelo_up, &
	971	zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
	972	zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
	973	)
[3999]	974
[4803]	975	ENDIF ! ok_poprecip
	976
[5408]	977	! End of precipitation processes after cloud formation
	978	! ----------------------------------------------------
[3999]	979
[5406]	980	!----------------------------------------------------------------------
[5408]	981	! P4> Calculation of cloud condensates amount seen by radiative scheme
[5406]	982	!----------------------------------------------------------------------
[4803]	983
[5413]	984	DO i=1,klon
[3999]	985
[5413]	986	IF (ok_poprecip) THEN
	987	IF (ok_radocond_snow) THEN
	988	radocond(i,k) = zoliq(i)
	989	zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
	990	ELSE
	991	radocond(i,k) = zoliq(i)
	992	zradoice(i) = zoliqi(i)
	993	ENDIF
	994	ELSE
	995	radocond(i,k) = zradocond(i)
	996	ENDIF
[4114]	997
[5413]	998	! calculate the percentage of ice in "radocond" so cloud+precip seen
	999	! by the radiation scheme
	1000	IF (radocond(i,k) .GT. 0.) THEN
	1001	radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
	1002	ENDIF
[4114]	1003	ENDDO
	1004
[3999]	1005	! ----------------------------------------------------------------
[5408]	1006	! P5> Wet scavenging
[3999]	1007	! ----------------------------------------------------------------
	1008
	1009	!Scavenging through nucleation in the layer
	1010
	1011	DO i = 1,klon
	1012
	1013	IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
	1014	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
	1015	ELSE
	1016	beta(i,k) = 0.
	1017	ENDIF
	1018
[4059]	1019	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
[3999]	1020
	1021	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
	1022
[4654]	1023	IF (temp(i,k) .GE. t_glace_min) THEN
[3999]	1024	zalpha_tr = a_tr_sca(3)
	1025	ELSE
	1026	zalpha_tr = a_tr_sca(4)
	1027	ENDIF
	1028
[5450]	1029	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
	1030	frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
[4226]	1031
[3999]	1032	! Nucleation with a factor of -1 instead of -0.5
[5450]	1033	zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
[3999]	1034
	1035	ENDIF
	1036
[4226]	1037	ENDDO
	1038
[3999]	1039	! Scavenging through impaction in the underlying layer
	1040
	1041	DO kk = k-1, 1, -1
	1042
	1043	DO i = 1, klon
	1044
	1045	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
	1046
[4654]	1047	IF (temp(i,kk) .GE. t_glace_min) THEN
[3999]	1048	zalpha_tr = a_tr_sca(1)
	1049	ELSE
	1050	zalpha_tr = a_tr_sca(2)
	1051	ENDIF
	1052
[5450]	1053	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
	1054	frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
[3999]	1055
	1056	ENDIF
	1057
	1058	ENDDO
	1059
	1060	ENDDO
[5406]	1061
	1062	!------------------------------------------------------------
[5408]	1063	! P6 > write diagnostics and outputs
[5406]	1064	!------------------------------------------------------------
	1065
	1066	!--AB Write diagnostics and tracers for ice supersaturation
	1067	IF ( ok_ice_supersat ) THEN
	1068	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
	1069	CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
[4226]	1070
[5406]	1071	DO i = 1, klon
	1072
[5545]	1073	cf_seri(i,k) = rneb(i,k)
	1074
[5406]	1075	IF ( zoliq(i) .LE. 0. ) THEN
	1076	!--If everything was precipitated, the remaining empty cloud is dissipated
	1077	!--and everything is transfered to the subsaturated clear sky region
[5545]	1078	!--NB. we do not change rneb, as it is a diagnostic only
	1079	cf_seri(i,k) = 0.
[5406]	1080	qvc(i) = 0.
	1081	ENDIF
	1082
	1083	IF ( .NOT. ok_unadjusted_clouds ) THEN
	1084	qvc(i) = zqs(i) * rneb(i,k)
	1085	ENDIF
	1086	IF ( zq(i) .GT. min_qParent ) THEN
	1087	rvc_seri(i,k) = qvc(i) / zq(i)
	1088	ELSE
[5456]	1089	rvc_seri(i,k) = min_ratio
[5406]	1090	ENDIF
	1091	!--The MIN barrier is NEEDED because of:
	1092	!-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
	1093	!-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
	1094	!--The MAX barrier is a safeguard that should not be activated
	1095	rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
[5452]	1096	IF ( ok_plane_contrail ) THEN
	1097	IF ( rneb(i,k) .GT. min_qParent ) THEN
[5456]	1098	rcont_seri(i,k) = contfra(i,k) / rneb(i,k)
[5452]	1099	ELSE
[5456]	1100	rcont_seri(i,k) = min_ratio
[5452]	1101	ENDIF
	1102	!--This barrier should never be activated
	1103	rcont_seri(i,k) = MIN(MAX(rcont_seri(i,k), 0.), 1.)
	1104	ENDIF
[5406]	1105
	1106	!--Diagnostics
	1107	gamma_cond(i,k) = gammasat(i)
	1108	subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
	1109	qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
[5412]	1110	qcld(i,k) = qvc(i) + zoliq(i)
[5406]	1111	ENDDO
	1112	ENDIF
	1113
[4803]	1114	! Outputs:
	1115	!-------------------------------
	1116	! Precipitation fluxes at layer interfaces
	1117	! + precipitation fractions +
	1118	! temperature and water species tendencies
	1119	DO i = 1, klon
	1120	psfl(i,k)=zifl(i)
	1121	prfl(i,k)=zrfl(i)
	1122	pfraclr(i,k)=znebprecipclr(i)
	1123	pfracld(i,k)=znebprecipcld(i)
	1124	d_q(i,k) = zq(i) - qt(i,k)
	1125	d_t(i,k) = zt(i) - temp(i,k)
[5412]	1126
	1127	IF (ok_bug_phase_lscp) THEN
	1128	d_ql(i,k) = (1-zfice(i))*zoliq(i)
	1129	d_qi(i,k) = zfice(i)*zoliq(i)
	1130	ELSE
	1131	d_ql(i,k) = zoliql(i)
	1132	d_qi(i,k) = zoliqi(i)
	1133	ENDIF
	1134
[4803]	1135	ENDDO
	1136
	1137
[5413]	1138	ENDDO ! loop on k from top to bottom
[4059]	1139
[3999]	1140
	1141	! Rain or snow at the surface (depending on the first layer temperature)
	1142	DO i = 1, klon
	1143	snow(i) = zifl(i)
	1144	rain(i) = zrfl(i)
[4392]	1145	! c_iso final output
[3999]	1146	ENDDO
	1147
	1148	IF (ncoreczq>0) THEN
	1149	WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
	1150	ENDIF
	1151
[4654]	1152
	1153	RETURN
	1154
[4380]	1155	END SUBROUTINE lscp
[3999]	1156	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1157
[4664]	1158	END MODULE lmdz_lscp

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

Download in other formats:

Original Format