Context Navigation

lmdz_lscp_condensation.F90 @ 5229

Last change on this file since 5229 was 5210, checked in by evignon, 8 weeks ago

changes in ice radius treatment in lmdz_lscp_condebsation + minor changes

Borella

File size: 46.9 KB

Rev	Line
[4951]	1	MODULE lmdz_lscp_condensation
	2	!----------------------------------------------------------------
	3	! Module for condensation of clouds routines
	4	! that are called in LSCP
	5
	6
	7	IMPLICIT NONE
	8
	9	CONTAINS
	10
	11	!**********************************************************************************
	12	SUBROUTINE condensation_lognormal( &
	13	klon, temp, qtot, qsat, gamma_cond, ratqs, &
	14	keepgoing, cldfra, qincld, qvc)
	15
	16	!----------------------------------------------------------------------
	17	! This subroutine calculates the formation of clouds, using a
	18	! statistical cloud scheme. It uses a generalised lognormal distribution
	19	! of total water in the gridbox
	20	! See Bony and Emanuel, 2001
	21	!----------------------------------------------------------------------
	22
	23	USE lmdz_lscp_ini, ONLY: eps
	24
	25	IMPLICIT NONE
	26
	27	!
	28	! Input
	29	!
	30	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
	31	!
	32	REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K]
	33	REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg]
	34	REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg]
	35	REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-]
	36	REAL, INTENT(IN) , DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-]
	37	LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed
	38	!
	39	! Output
	40	! NB. those are in INOUT because of the convergence loop on temperature
	41	! (in some cases, the values are not re-computed) but the values
	42	! are never used explicitely
	43	!
	44	REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-]
	45	REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg]
	46	REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg]
	47	!
	48	! Local
	49	!
	50	INTEGER :: i
	51	REAL :: pdf_std, pdf_k, pdf_delta
	52	REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
	53	!
	54	!--Loop on klon
	55	!
	56	DO i = 1, klon
	57
	58	!--If a new calculation of the condensation is needed,
	59	!--i.e., temperature has not yet converged (or the cloud is
	60	!--formed elsewhere)
	61	IF (keepgoing(i)) THEN
	62
	63	pdf_std = ratqs(i) * qtot(i)
	64	pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
	65	pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
	66	pdf_a = pdf_delta / ( pdf_k * SQRT(2.) )
	67	pdf_b = pdf_k / (2. * SQRT(2.))
	68	pdf_e1 = pdf_a - pdf_b
	69	pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
	70	pdf_e1 = 1. - ERF(pdf_e1)
	71	pdf_e2 = pdf_a + pdf_b
	72	pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
	73	pdf_e2 = 1. - ERF(pdf_e2)
	74
	75
	76	IF ( pdf_e1 .LT. eps ) THEN
	77	cldfra(i) = 0.
	78	qincld(i) = qsat(i)
	79	!--AB grid-mean vapor in the cloud - we assume saturation adjustment
	80	qvc(i) = 0.
	81	ELSE
	82	cldfra(i) = 0.5 * pdf_e1
	83	qincld(i) = qtot(i) * pdf_e2 / pdf_e1
	84	!--AB grid-mean vapor in the cloud - we assume saturation adjustment
	85	qvc(i) = qsat(i) * cldfra(i)
	86	ENDIF
	87
	88	ENDIF ! end keepgoing
	89
	90	ENDDO ! end loop on i
	91	END SUBROUTINE condensation_lognormal
	92	!**********************************************************************************
	93
	94	!**********************************************************************************
	95	SUBROUTINE condensation_ice_supersat( &
	96	klon, dtime, missing_val, pplay, paprsdn, paprsup, &
	97	cf_seri, rvc_seri, ratio_qi_qtot, shear, pbl_eps, cell_area, &
	98	temp, qtot, qsat, gamma_cond, ratqs, keepgoing, &
	99	cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
	100	dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
	101	Tcontr, qcontr, qcontr2, fcontrN, fcontrP, flight_dist, flight_h2o, &
	102	dcf_avi, dqi_avi, dqvc_avi)
	103
	104	!----------------------------------------------------------------------
	105	! This subroutine calculates the formation, evolution and dissipation
	106	! of clouds, using a process-oriented treatment of the cloud properties
	107	! (cloud fraction, vapor in the cloud, condensed water in the cloud).
	108	! It allows for ice supersaturation in cold regions, in clear sky.
	109	! If ok_unadjusted_clouds, it also allows for sub- and supersaturated
	110	! cloud water vapors.
	111	! It also allows for the formation and evolution of condensation trails
	112	! (contrails) from aviation.
	113	! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
	114	! April 2024
	115	!----------------------------------------------------------------------
	116
	117	USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
	118	USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
	119	USE lmdz_lscp_ini, ONLY: eps, temp_nowater, ok_weibull_warm_clouds
	120	USE lmdz_lscp_ini, ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf
	121	USE lmdz_lscp_ini, ONLY: lunout
	122
[5165]	123	USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, &
	124	mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
[4951]	125	rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, &
[5210]	126	coef_mixing_lscp, coef_shear_lscp, &
[4951]	127	chi_mixing_lscp, rho_ice
	128
	129	IMPLICIT NONE
	130
	131	!
	132	! Input
	133	!
	134	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
	135	REAL, INTENT(IN) :: dtime ! time step [s]
	136	REAL, INTENT(IN) :: missing_val ! missing value for output
	137	!
	138	REAL, INTENT(IN) , DIMENSION(klon) :: pplay ! layer pressure [Pa]
	139	REAL, INTENT(IN) , DIMENSION(klon) :: paprsdn ! pressure at the lower interface [Pa]
	140	REAL, INTENT(IN) , DIMENSION(klon) :: paprsup ! pressure at the upper interface [Pa]
	141	REAL, INTENT(IN) , DIMENSION(klon) :: cf_seri ! cloud fraction [-]
	142	REAL, INTENT(IN) , DIMENSION(klon) :: rvc_seri ! gridbox-mean water vapor in cloud [kg/kg]
	143	REAL, INTENT(IN) , DIMENSION(klon) :: ratio_qi_qtot ! specific ice water content to total specific water ratio [-]
	144	REAL, INTENT(IN) , DIMENSION(klon) :: shear ! vertical shear [s-1]
	145	REAL, INTENT(IN) , DIMENSION(klon) :: pbl_eps !
	146	REAL, INTENT(IN) , DIMENSION(klon) :: cell_area !
	147	REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K]
	148	REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg]
	149	REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg]
	150	REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-]
	151	REAL, INTENT(INOUT), DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-]
	152	LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed
	153	!
	154	! Input for aviation
	155	!
	156	REAL, INTENT(IN), DIMENSION(klon) :: flight_dist !
	157	REAL, INTENT(IN), DIMENSION(klon) :: flight_h2o !
	158	!
	159	! Output
	160	! NB. cldfra and qincld should be outputed as cf_seri and qi_seri,
	161	! or as tendencies (maybe in the future)
	162	! NB. those are in INOUT because of the convergence loop on temperature
	163	! (in some cases, the values are not re-computed) but the values
	164	! are never used explicitely
	165	!
	166	REAL, INTENT(INOUT), DIMENSION(klon) :: cldfra ! cloud fraction [-]
	167	REAL, INTENT(INOUT), DIMENSION(klon) :: qincld ! cloud-mean in-cloud total specific water [kg/kg]
	168	REAL, INTENT(INOUT), DIMENSION(klon) :: qvc ! gridbox-mean vapor in the cloud [kg/kg]
	169	REAL, INTENT(INOUT), DIMENSION(klon) :: issrfra ! ISSR fraction [-]
	170	REAL, INTENT(INOUT), DIMENSION(klon) :: qissr ! gridbox-mean ISSR specific water [kg/kg]
	171	!
	172	! Diagnostics for condensation and ice supersaturation
	173	! NB. idem for the INOUT
	174	!
	175	REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_sub ! cloud fraction tendency because of sublimation [s-1]
	176	REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_con ! cloud fraction tendency because of condensation [s-1]
	177	REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_mix ! cloud fraction tendency because of cloud mixing [s-1]
	178	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_adj ! specific ice content tendency because of temperature adjustment [kg/kg/s]
	179	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_sub ! specific ice content tendency because of sublimation [kg/kg/s]
	180	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_con ! specific ice content tendency because of condensation [kg/kg/s]
	181	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_mix ! specific ice content tendency because of cloud mixing [kg/kg/s]
	182	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
	183	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
	184	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
	185	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
	186	!
	187	! Diagnostics for aviation
	188	! NB. idem for the INOUT
	189	!
	190	REAL, INTENT(INOUT), DIMENSION(klon) :: Tcontr ! critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) [K]
	191	REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr !
	192	REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr2 !
	193	REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrN !
	194	REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrP !
	195	REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_avi ! cloud fraction tendency because of aviation [s-1]
	196	REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_avi ! specific ice content tendency because of aviation [kg/kg/s]
	197	REAL, INTENT(INOUT), DIMENSION(klon) :: dqvc_avi ! specific cloud water vapor tendency because of aviation [kg/kg/s]
	198	!
	199	! Local
	200	!
	201	INTEGER :: i
	202	LOGICAL :: ok_warm_cloud
	203	REAL, DIMENSION(klon) :: qcld, qzero
	204	!
	205	! for lognormal
	206	REAL :: pdf_std, pdf_k, pdf_delta
	207	REAL :: pdf_a, pdf_b, pdf_e1, pdf_e2
	208	!
	209	! for unadjusted clouds
	210	REAL :: qvapincld, qvapincld_new
	211	REAL :: qiceincld, qice_ratio
	212	REAL :: pres_sat, kappa
	213	REAL :: air_thermal_conduct, water_vapor_diff
[5210]	214	REAL :: iwc
	215	REAL :: iwc_log_inf100, iwc_log_sup100
	216	REAL :: iwc_inf100, alpha_inf100, coef_inf100
	217	REAL :: mu_sup100, sigma_sup100, coef_sup100
	218	REAL :: Dm_ice, rm_ice
[4951]	219	!
	220	! for sublimation
	221	REAL :: pdf_alpha
	222	REAL :: dqt_sub
	223	!
	224	! for condensation
	225	REAL, DIMENSION(klon) :: qsatl, dqsatl
	226	REAL :: clrfra, qclr, sl_clr, rhl_clr
	227	REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc
	228	REAL :: pdf_x, pdf_y, pdf_T
	229	REAL :: pdf_e3, pdf_e4
	230	REAL :: cf_cond, qt_cond, dqt_con
	231	!
	232	! for mixing
	233	REAL, DIMENSION(klon) :: subfra, qsub
	234	REAL :: dqt_mix_sub, dqt_mix_issr
	235	REAL :: dcf_mix_sub, dcf_mix_issr
	236	REAL :: dqvc_mix_sub, dqvc_mix_issr
	237	REAL :: dqt_mix
	238	REAL :: a_mix, bovera, N_cld_mix, L_mix
	239	REAL :: envfra_mix, cldfra_mix
	240	REAL :: L_shear, shear_fra
	241	REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix
	242	!
	243	! for cell properties
	244	REAL :: rho, rhodz, dz
	245	!REAL :: V_cell, M_cell
	246	!
	247	! for aviation and cell properties
	248	!REAL :: dqt_avi
	249	!REAL :: contrail_fra
	250	!
	251	!
	252	!--more local variables for diagnostics
	253	!--imported from YOMCST.h
	254	!--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1)
	255	!--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air
	256	!--values from Schumann, Meteorol Zeitschrift, 1996
	257	!--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
	258	!--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
	259	!REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1)
	260	!REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1)
	261	!REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft
	262	!--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
	263	!--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
	264	!REAL :: Gcontr
	265	!--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2)
	266	!--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1)
	267	!REAL :: qsatliqcontr
	268
	269
	270	!-----------------------------------------------
	271	! Initialisations
	272	!-----------------------------------------------
	273
	274	! Ajout des émissions de H2O dues à l'aviation
	275	! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
	276	! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O )
	277	! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
	278	! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
	279	! flight_h2O is in kg H2O / s / cell
	280	!
	281	!IF (ok_plane_h2o) THEN
	282	! q = ( M_cellq + flight_h2o(i,k)dtime(1.-q) ) / (M_cell + flight_h2o(i,k)dtime*(1.-q) )
	283	!ENDIF
	284
	285
	286	qzero(:) = 0.
	287
	288	!--Calculation of qsat w.r.t. liquid
	289	CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl)
	290
	291	!
	292	!--Loop on klon
	293	!
	294	DO i = 1, klon
	295
	296	!--If a new calculation of the condensation is needed,
	297	!--i.e., temperature has not yet converged (or the cloud is
	298	!--formed elsewhere)
	299	IF (keepgoing(i)) THEN
	300
	301	!--If the temperature is higher than the threshold below which
	302	!--there is no liquid in the gridbox, we activate the usual scheme
	303	!--(generalised lognormal from Bony and Emanuel 2001)
	304	!--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
	305	!--all clouds, and the lognormal scheme is not activated
	306	IF ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) THEN
	307
	308	pdf_std = ratqs(i) * qtot(i)
	309	pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2. ) )
	310	pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
	311	pdf_a = pdf_delta / ( pdf_k * SQRT(2.) )
	312	pdf_b = pdf_k / (2. * SQRT(2.))
	313	pdf_e1 = pdf_a - pdf_b
	314	pdf_e1 = SIGN( MIN(ABS(pdf_e1), 5.), pdf_e1 )
	315	pdf_e1 = 1. - ERF(pdf_e1)
	316	pdf_e2 = pdf_a + pdf_b
	317	pdf_e2 = SIGN( MIN(ABS(pdf_e2), 5.), pdf_e2 )
	318	pdf_e2 = 1. - ERF(pdf_e2)
	319
	320
	321	IF ( pdf_e1 .LT. eps ) THEN
	322	cldfra(i) = 0.
	323	qincld(i) = qsat(i)
	324	qvc(i) = 0.
	325	ELSE
	326	cldfra(i) = 0.5 * pdf_e1
	327	qincld(i) = qtot(i) * pdf_e2 / pdf_e1
	328	qvc(i) = qsat(i) * cldfra(i)
	329	ENDIF
	330
	331	!--If the temperature is lower than temp_nowater, we use the new
	332	!--condensation scheme that allows for ice supersaturation
	333	ELSE
	334
	335	!--Initialisation
	336	IF ( temp(i) .GT. temp_nowater ) THEN
	337	!--If the air mass is warm (liquid water can exist),
	338	!--all the memory is lost and the scheme becomes statistical,
	339	!--i.e., the sublimation and mixing processes are deactivated,
	340	!--and the condensation process is slightly adapted
	341	!--This can happen only if ok_weibull_warm_clouds = .TRUE.
	342	cldfra(i) = 0.
	343	qvc(i) = 0.
	344	qcld(i) = 0.
	345	ok_warm_cloud = .TRUE.
	346	ELSE
	347	!--The following barriers ensure that the traced cloud properties
	348	!--are consistent. In some rare cases, i.e. the cloud water vapor
	349	!--can be greater than the total water in the gridbox
	350	cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
	351	qcld(i) = MAX(0., MIN(qtot(i), ( ratio_qi_qtot(i) + rvc_seri(i) ) * qtot(i)))
	352	qvc(i) = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
	353	ok_warm_cloud = .FALSE.
	354	ENDIF
	355
	356	dcf_sub(i) = 0.
	357	dqi_sub(i) = 0.
	358	dqvc_sub(i) = 0.
	359	dqi_adj(i) = 0.
	360	dqvc_adj(i) = 0.
	361	dcf_con(i) = 0.
	362	dqi_con(i) = 0.
	363	dqvc_con(i) = 0.
	364	dcf_mix(i) = 0.
	365	dqi_mix(i) = 0.
	366	dqvc_mix(i) = 0.
	367
	368	issrfra(i) = 0.
	369	qissr(i) = 0.
	370	subfra(i) = 0.
	371	qsub(i) = 0.
	372
	373	!--Initialisation of the cell properties
	374	!--Dry density [kg/m3]
	375	rho = pplay(i) / temp(i) / RD
	376	!--Dry air mass [kg/m2]
	377	rhodz = ( paprsdn(i) - paprsup(i) ) / RG
	378	!--Cell thickness [m]
	379	dz = rhodz / rho
	380	!--Cell volume [m3]
	381	!V_cell = dz * cell_area(i)
	382	!--Cell dry air mass [kg]
	383	!M_cell = rhodz * cell_area(i)
	384
	385
	386	IF ( ok_unadjusted_clouds ) THEN
	387	!--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
	388	!--hypothesis is lost, and the vapor in the cloud is purely prognostic.
	389	!
	390	!--The deposition equation is
[5041]	391	!-- dmi/dt = alpha4piCSvi / ( R_vT/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
	392	!--from Lohmann et al. (2016), where
	393	!--alpha is the deposition coefficient [-]
[4951]	394	!--mi is the mass of one ice crystal [kg]
	395	!--C is the capacitance of an ice crystal [m]
	396	!--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
[5041]	397	!--R_v is the specific gas constant for humid air [J/kg/K]
[4951]	398	!--T is the temperature [K]
	399	!--esi is the saturation pressure w.r.t. ice [Pa]
	400	!--Dv is the diffusivity of water vapor [m2/s]
	401	!--Ls is the specific latent heat of sublimation [J/kg/K]
	402	!--ka is the thermal conductivity of dry air [J/m/s/K]
	403	!
[5165]	404	!--alpha is a coefficient to take into account the fact that during deposition, a water
	405	!--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
	406	!--coherent (with the same structure). It has no impact for sublimation.
	407	!--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
	408	!--during deposition, and alpha = 1. during sublimation.
[5041]	409	!--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
[5210]	410	!-- C = capa_cond_cirrus * rm_ice
[5041]	411	!
[4951]	412	!--We have qice = Nice * mi, where Nice is the ice crystal
	413	!--number concentration per kg of moist air
	414	!--HYPOTHESIS 1: the ice crystals are spherical, therefore
[5210]	415	!-- mi = 4/3 * pi * rm_ice*3 rho_ice
[4951]	416	!--HYPOTHESIS 2: the ice crystals are monodisperse with the
[5210]	417	!--initial radius rm_ice_0.
[4951]	418	!--NB. this is notably different than the assumption
	419	!--of a distributed qice in the cloud made in the sublimation process;
	420	!--should it be consistent?
	421	!
	422	!--As the deposition process does not create new ice crystals,
[5210]	423	!--and because we assume a same rm_ice value for all crystals
[4951]	424	!--therefore the sublimation process does not destroy ice crystals
	425	!--(or, in a limit case, it destroys all ice crystals), then
	426	!--Nice is a constant during the sublimation/deposition process.
[5210]	427	!-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
[4951]	428	!
	429	!--The deposition equation then reads:
[5210]	430	!-- dqi/dt = alpha4picapa_cond_cirrusrm_ice(qvc-qsat)/qsat / ( R_vT/esi/Dv + Ls/ka/T (Ls/R_v/T - 1) ) * Nice
	431	!-- dqi/dt = alpha4picapa_cond_cirrus* (qi / qi_0)*(1/3) rm_ice_0*(qvc-qsat)/qsat &
[4951]	432	!-- / ( R_vT/esi/Dv + Ls/ka/T (Ls*R_v/T - 1) ) &
[5210]	433	!-- * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
[4951]	434	!-- dqi/dt = qi*(1/3) (qvc - qsat) * qi_0**(2/3) &
[5210]	435	!-- alpha/qsatcapa_cond_cirrus/ (R_vT/esi/Dv + Ls/ka/T(LsR_v/T - 1)) / ( 1/3 rm_ice_0*2 rho_ice )
[4951]	436	!--and we have
[5210]	437	!-- dqvc/dt = - qi*(1/3) (qvc - qsat) / kappa * alpha * qi_0(2/3) / rm_ice_02
	438	!-- dqi/dt = qi*(1/3) (qvc - qsat) / kappa * alpha * qi_0(2/3) / rm_ice_02
[5165]	439	!--where kappa = 1/3rho_ice/capa_cond_cirrusqsat(R_vT/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
[4951]	440	!
	441	!--This system of equations can be resolved with an exact
	442	!--explicit numerical integration, having one variable resolved
	443	!--explicitly, the other exactly. The exactly resolved variable is
	444	!--the one decreasing, so it is qvc if the process is
	445	!--condensation, qi if it is sublimation.
	446	!
	447	!--kappa is computed as an initialisation constant, as it depends only
	448	!--on temperature and other pre-computed values
	449	pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i)
	450	!--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
	451	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
	452	!--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
	453	water_vapor_diff = 0.211 * ( temp(i) / RTT )*1.94 ( 101325. / pplay(i) ) * 1.e-4
[5165]	454	kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat(i) &
[4951]	455	* ( RV * temp(i) / water_vapor_diff / pres_sat &
	456	+ RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
[5041]	457	!--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
[4951]	458	ENDIF
	459
	460
	461	!-------------------------------------------------------------------
	462	!-- SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD --
	463	!-------------------------------------------------------------------
	464
	465	!--If there is a cloud
	466	IF ( cldfra(i) .GT. eps ) THEN
	467
	468	qvapincld = qvc(i) / cldfra(i)
	469	qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
	470
	471	!--If the ice water content is too low, the cloud is purely sublimated
[4974]	472	!--Most probably, we advected a cloud with no ice water content (possible
	473	!--if the entire cloud precipited for example)
[4951]	474	IF ( qiceincld .LT. eps ) THEN
	475	dcf_sub(i) = - cldfra(i)
	476	dqvc_sub(i) = - qvc(i)
	477	dqi_sub(i) = - ( qcld(i) - qvc(i) )
	478
	479	cldfra(i) = 0.
	480	qcld(i) = 0.
	481	qvc(i) = 0.
	482
	483	!--Else, the cloud is adjusted and sublimated
	484	ELSE
	485
[4974]	486	!--The vapor in cloud cannot be higher than the
	487	!--condensation threshold
	488	qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
	489	qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
	490
[4951]	491	IF ( ok_unadjusted_clouds ) THEN
	492	!--Here, the initial vapor in the cloud is qvapincld, and we compute
	493	!--the new vapor qvapincld_new
	494
[5210]	495	!--rm_ice formula from McFarquhar and Heymsfield (1997)
	496	iwc = qiceincld * rho * 1e3
	497	iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
	498	iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
	499	iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
[4951]	500
[5210]	501	alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
	502	coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120.
	503
	504	mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) &
	505	+ ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100
	506	sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) &
	507	+ ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100
	508	coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3 * mu_sup100 + 4.5 * sigma_sup100**2. )
	509
	510	Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) &
	511	* coef_sup100 ) / ( coef_inf100 + coef_sup100 )
	512	rm_ice = Dm_ice / 2. * 1.E-6
	513
[4951]	514	IF ( qvapincld .GE. qsat(i) ) THEN
	515	!--If the cloud is initially supersaturated
	516	!--Exact explicit integration (qvc exact, qice explicit)
	517	qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) &
[5210]	518	* EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2. )
[4951]	519	ELSE
	520	!--If the cloud is initially subsaturated
	521	!--Exact explicit integration (qice exact, qvc explicit)
	522	!--The barrier is set so that the resulting vapor in cloud
	523	!--cannot be greater than qsat
	524	!--qice_ratio is the ratio between the new ice content and
	525	!--the old one, it is comprised between 0 and 1
[5210]	526	qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice*2. dtime * ( qsat(i) - qvapincld ) )
[4951]	527
	528	IF ( qice_ratio .LT. 0. ) THEN
	529	!--If all the ice has been sublimated, we sublimate
	530	!--completely the cloud and do not activate the sublimation
	531	!--process
	532	!--Tendencies and diagnostics
	533	dcf_sub(i) = - cldfra(i)
	534	dqvc_sub(i) = - qvc(i)
	535	dqi_sub(i) = - ( qcld(i) - qvc(i) )
	536
	537	cldfra(i) = 0.
	538	qcld(i) = 0.
	539	qvc(i) = 0.
	540
	541	!--The new vapor in cloud is set to 0 so that the
	542	!--sublimation process does not activate
	543	qvapincld_new = 0.
	544	ELSE
	545	!--Else, the sublimation process is activated with the
	546	!--diagnosed new cloud water vapor
	547	!--The new vapor in the cloud is increased with the
	548	!--sublimated ice
	549	qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**(3./2.) )
	550	!--The new vapor in the cloud cannot be greater than qsat
	551	qvapincld_new = MIN(qvapincld_new, qsat(i))
	552	ENDIF ! qice_ratio .LT. 0.
	553	ENDIF ! qvapincld .GT. qsat(i)
	554	ELSE
	555	!--We keep the saturation adjustment hypothesis, and the vapor in the
	556	!--cloud is set equal to the saturation vapor
	557	qvapincld_new = qsat(i)
	558	ENDIF ! ok_unadjusted_clouds
[5210]	559
	560	!--Adjustment of the IWC to the new vapor in cloud
	561	!--(this can be either positive or negative)
	562	dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
	563	dqi_adj(i) = - dqvc_adj(i)
[4951]	564
[5210]	565	!--Add tendencies
	566	!--The vapor in the cloud is updated, but not qcld as it is constant
	567	!--through this process, as well as cldfra which is unmodified
	568	qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
	569
	570
	571	!------------------------------------
	572	!-- DISSIPATION OF THE CLOUD --
	573	!------------------------------------
	574
[4951]	575	!--If the vapor in cloud is below vapor needed for the cloud to survive
	576	IF ( qvapincld .LT. qvapincld_new ) THEN
	577	!--Sublimation of the subsaturated cloud
	578	!--iflag_cloud_sublim_pdf selects the PDF of the ice water content
	579	!--to use.
	580	!--iflag = 1 --> uniform distribution
	581	!--iflag = 2 --> exponential distribution
	582	!--iflag = 3 --> gamma distribution (Karcher et al 2018)
	583
	584	IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
	585	!--Uniform distribution starting at qvapincld
	586	pdf_e1 = 1. / ( 2. * qiceincld )
	587
	588	dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
	589	dqt_sub = - cldfra(i) * ( qvapincld_new2. - qvapincld2. ) &
	590	* pdf_e1 / 2.
	591
	592	ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
	593	!--Exponential distribution starting at qvapincld
	594	pdf_alpha = 1. / qiceincld
	595	pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
	596
	597	dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
	598	dqt_sub = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
	599	+ qvapincld - qvapincld_new * pdf_e1 )
	600
	601	ELSEIF ( iflag_cloud_sublim_pdf .GE. 3 ) THEN
	602	!--Gamma distribution starting at qvapincld
	603	pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
	604	pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
	605	pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
	606	pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
	607
	608	dcf_sub(i) = - cldfra(i) * pdf_e1
	609	dqt_sub = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
	610	ENDIF
	611
	612	!--Tendencies and diagnostics
[5210]	613	dqvc_sub(i) = dqt_sub
[4951]	614
	615	!--Add tendencies
	616	cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
	617	qcld(i) = MAX(0., qcld(i) + dqt_sub)
	618	qvc(i) = MAX(0., qvc(i) + dqvc_sub(i))
	619
	620	ENDIF ! qvapincld .LT. qvapincld_new
	621
	622	ENDIF ! qiceincld .LT. eps
	623	ENDIF ! cldfra(i) .GT. eps
	624
	625
	626	!--------------------------------------------------------------------------
	627	!-- CONDENSATION AND DIAGNOTICS OF SUB- AND SUPERSATURATED REGIONS --
	628	!--------------------------------------------------------------------------
	629	!--This section relies on a distribution of water in the clear-sky region of
	630	!--the mesh.
	631
	632	!--If there is a clear-sky region
	633	IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
	634
	635	!--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
	636	qclr = qtot(i) - qcld(i)
	637
	638	!--New PDF
	639	rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
	640	rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
	641
	642	!--Calculation of the properties of the PDF
	643	!--Parameterization from IAGOS observations
	644	!--pdf_e1 and pdf_e2 will be reused below
	645
[5161]	646	!--Coefficient for standard deviation:
	647	!-- tuning coef * (clear sky area*0.25) (function of temperature)
	648	pdf_e1 = beta_pdf_lscp &
	649	* ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
[4951]	650	* MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
	651	IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
	652	pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
	653	ELSE
	654	pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
	655	ENDIF
	656	pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
	657	pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
	658
	659	pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
	660	pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
	661	pdf_loc = rhl_clr - pdf_scale * pdf_e2
	662
	663	!--Diagnostics of ratqs
	664	ratqs(i) = pdf_std / ( qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. )
	665
	666	!--Calculation of the newly condensed water and fraction (pronostic)
	667	!--Integration of the clear sky PDF between gamma_cond*qsat and +inf
	668	!--NB. the calculated values are clear-sky averaged
	669
	670	pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
	671	pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
	672	pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
	673	pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
	674	cf_cond = EXP( - pdf_y )
	675	qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
	676
	677	IF ( ok_warm_cloud ) THEN
[5210]	678	!--If the statistical scheme is activated, the calculated increase is equal
	679	!--to the cloud fraction, we assume saturation adjustment, and the
	680	!--condensation process stops
[4951]	681	cldfra(i) = cf_cond
	682	qcld(i) = qt_cond
	683	qvc(i) = cldfra(i) * qsat(i)
	684
[5210]	685	ELSEIF ( cf_cond .GT. eps ) THEN
[4951]	686
[5210]	687	dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
	688	dqt_con = ( 1. - cldfra(i) ) * qt_cond
[4951]	689
	690	!--Barriers
	691	dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
	692	dqt_con = MIN(dqt_con, qclr)
	693
	694
	695	IF ( ok_unadjusted_clouds ) THEN
	696	!--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute
	697	!--the new vapor qvapincld. The timestep is divided by two because we do not
	698	!--know when the condensation occurs
	699	qvapincld = gamma_cond(i) * qsat(i)
	700	qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
[5210]	701
	702	!--rm_ice formula from McFarquhar and Heymsfield (1997)
	703	iwc = qiceincld * rho * 1e3
	704	iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
	705	iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
	706	iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
	707
	708	alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
	709	coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120.
	710
	711	mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) &
	712	+ ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100
	713	sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) &
	714	+ ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100
	715	coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2. )
	716
	717	Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) &
	718	* coef_sup100 ) / ( coef_inf100 + coef_sup100 )
	719	rm_ice = Dm_ice / 2. * 1.E-6
[4951]	720	!--As qvapincld is necessarily greater than qsat, we only
	721	!--use the exact explicit formulation
	722	!--Exact explicit version
	723	qvapincld = qsat(i) + ( qvapincld - qsat(i) ) &
[5210]	724	* EXP( - depo_coef_cirrus * dtime / 2. * qiceincld / kappa / rm_ice**2. )
[4951]	725	ELSE
	726	!--We keep the saturation adjustment hypothesis, and the vapor in the
	727	!--newly formed cloud is set equal to the saturation vapor.
	728	qvapincld = qsat(i)
	729	ENDIF
	730
	731	!--Tendency on cloud vapor and diagnostic
	732	dqvc_con(i) = qvapincld * dcf_con(i)
	733	dqi_con(i) = dqt_con - dqvc_con(i)
	734
	735	!--Note that the tendencies are NOT added because they are
	736	!--added after the mixing process. In the following, the gridbox fraction is
	737	!-- 1. - dcf_con(i), and the total water in the gridbox is
	738	!-- qtot(i) - dqi_con(i) - dqvc_con(i)
	739
[5210]	740	ENDIF ! ok_warm_cloud, cf_cond .GT. eps
[4951]	741	ENDIF ! ( 1. - cldfra(i) ) .GT. eps
	742
	743	!--If there is still clear sky, we diagnose the ISSR
	744	!--We recalculte the PDF properties (after the condensation process)
	745	IF ( ( ( 1. - dcf_con(i) - cldfra(i) ) .GT. eps ) .AND. .NOT. ok_warm_cloud ) THEN
	746	!--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
	747	qclr = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i)
	748
	749	!--New PDF
	750	rhl_clr = qclr / ( 1. - dcf_con(i) - cldfra(i) ) / qsatl(i) * 100.
	751	rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)
	752
	753	!--Calculation of the properties of the PDF
	754	!--Parameterization from IAGOS observations
	755	!--pdf_e1 and pdf_e2 will be reused below
	756
[5161]	757	!--Coefficient for standard deviation:
	758	!-- tuning coef * (clear sky area*0.25) (function of temperature)
	759	pdf_e1 = beta_pdf_lscp &
	760	* ( ( 1. - dcf_con(i) - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
[4951]	761	* MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
	762	IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN
	763	pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp
	764	ELSE
	765	pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp
	766	ENDIF
	767	pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
	768	pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3
	769
	770	pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
	771	pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ))
	772	pdf_loc = rhl_clr - pdf_scale * pdf_e2
	773
	774	!--We then calculate the part that is greater than qsat
	775	!--and consider it supersaturated
	776
	777	pdf_x = qsat(i) / qsatl(i) * 100.
	778	pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
	779	pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
	780	pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
	781	issrfra(i) = EXP( - pdf_y ) * ( 1. - dcf_con(i) - cldfra(i) )
	782	qissr(i) = ( pdf_e3 * ( 1. - dcf_con(i) - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
	783	ENDIF
	784
	785	!--Calculation of the subsaturated clear sky fraction and water
	786	subfra(i) = 1. - dcf_con(i) - cldfra(i) - issrfra(i)
	787	qsub(i) = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i) - qissr(i)
	788
	789
	790	!--------------------------------------
	791	!-- CLOUD MIXING --
	792	!--------------------------------------
	793	!--This process mixes the cloud with its surroundings: the subsaturated clear sky,
	794	!--and the supersaturated clear sky. It is activated if the cloud is big enough,
	795	!--but does not cover the entire mesh.
	796	!
	797	IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) &
	798	.AND. .NOT. ok_warm_cloud ) THEN
	799
	800	!--Initialisation
	801	dcf_mix_sub = 0.
	802	dqt_mix_sub = 0.
	803	dqvc_mix_sub = 0.
	804	dcf_mix_issr = 0.
	805	dqt_mix_issr = 0.
	806	dqvc_mix_issr = 0.
	807
	808
	809	!-- PART 1 - TURBULENT DIFFUSION
	810
	811	!--Clouds within the mesh are assumed to be ellipses. The length of the
	812	!--semi-major axis is a and the length of the semi-minor axis is b.
	813	!--N_cld_mix is the number of clouds within the mesh, and
	814	!--clouds_perim is the total perimeter of the clouds within the mesh,
	815	!--not considering interfaces with other meshes (only the interfaces with clear
	816	!--sky are taken into account).
	817	!--
	818	!--The area of each cloud is A = a * b * RPI,
	819	!--and the perimeter of each cloud is
	820	!-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
	821	!--
	822	!--With cell_area the area of the cell, we have:
	823	!-- cldfra = A * N_cld_mix / cell_area
	824	!-- clouds_perim = P * N_cld_mix
	825	!--
	826	!--We assume that the ratio between b and a is a function of
	827	!--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
	828	!--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
	829	!--are spherical.
	830	!-- b / a = bovera = MAX(0.1, cldfra)
	831	bovera = MAX(0.1, cldfra(i))
	832	!--The clouds perimeter is imposed using the formula from Morcrette 2012,
	833	!--based on observations.
	834	!-- clouds_perim_normalized = alpha * cldfra * ( 1. - cldfra ) = clouds_perim / ( norm_constant * cell_area )
	835	!--
	836	!--With all this, we have
	837	!-- a = SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) )
	838	!-- P = RPI * a * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
	839	!--and therefore, using the perimeter
	840	!-- alpha * cldfra * ( 1. - cldfra ) * norm_constant * cell_area
	841	!-- = N_cld_mix * RPI &
	842	!-- * SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) &
	843	!-- * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
	844	!--and finally
	845	N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - dcf_con(i) - cldfra(i) )**2. &
	846	* cell_area(i) * ( 1. - dcf_con(i) ) * bovera / RPI &
	847	/ ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2.
	848	!--where coef_mix_lscp = ( alpha * norm_constant )**2.
	849	!--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer
	850	!--In particular, it is 0 if cldfra = 1
	851	a_mix = SQRT( cell_area(i) * ( 1. - dcf_con(i) ) * cldfra(i) / bovera / N_cld_mix / RPI )
	852
	853	!--The time required for turbulent diffusion to homogenize a region of size
	854	!--L_mix is defined as (L_mix2./tke_dissip)(1./3.) (Pope, 2000; Field et al., 2014)
	855	!--We compute L_mix and assume that the cloud is mixed over this length
	856	L_mix = SQRT( dtime*3. pbl_eps(i) )
	857	!--The mixing length cannot be greater than the semi-minor axis. In this case,
	858	!--the entire cloud is mixed.
	859	L_mix = MIN(L_mix, a_mix * bovera)
	860
	861	!--The fraction of clear sky mixed is
	862	!-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
	863	envfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
	864	* ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2. )
	865	!--The fraction of cloudy sky mixed is
	866	!-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
	867	cldfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) ) &
	868	* ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2. )
	869
	870
	871	!-- PART 2 - SHEARING
	872
	873	!--The clouds are then sheared. We keep the shape and number
	874	!--assumptions from before. The clouds are sheared along their
	875	!--semi-major axis (a_mix), on the entire cell heigh dz.
	876	!--The increase in size is
	877	L_shear = coef_shear_lscp * shear(i) * dz * dtime
	878	!--therefore, the total increase in fraction is
[5210]	879	!-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
	880	shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix &
[4951]	881	/ cell_area(i) / ( 1. - dcf_con(i) )
	882	!--and the environment and cloud mixed fractions are the same,
	883	!--which we add to the previous calculated mixed fractions.
	884	!--We therefore assume that the sheared clouds and the turbulent
	885	!--mixed clouds are different.
	886	envfra_mix = envfra_mix + shear_fra
	887	cldfra_mix = cldfra_mix + shear_fra
	888
	889
	890	!-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
	891
	892	!--The environment fraction is allocated to subsaturated sky or supersaturated sky,
	893	!--according to the factor sigma_mix. This is computed as the ratio of the
	894	!--subsaturated sky fraction to the environment fraction, corrected by a factor
	895	!--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
	896	!--supersaturated sky is favoured. Physically, this means that it is more likely
	897	!--to have supersaturated sky around the cloud than subsaturated sky.
	898	sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
	899	subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
	900	issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
	901	cldfra_mix = MIN( cldfra_mix, cldfra(i) )
	902
	903	!--First, we mix the subsaturated sky (subfra_mix) and the cloud close
	904	!--to this fraction (sigma_mix * cldfra_mix).
	905	IF ( subfra(i) .GT. eps ) THEN
	906
	907	IF ( ok_unadjusted_clouds ) THEN
	908	!--The subsaturated air is simply added to the cloud,
	909	!--with the corresponding cloud fraction
	910	!--If the cloud is too subsaturated, the sublimation process
	911	!--activated in the following timestep will reduce the cloud
	912	!--fraction
	913	dcf_mix_sub = subfra_mix
	914	dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
	915	dqvc_mix_sub = dqt_mix_sub
	916
	917	ELSE
	918	!--We compute the total humidity in the mixed air, which
	919	!--can be either sub- or supersaturated.
	920	qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
	921	+ qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
	922	/ ( subfra_mix + cldfra_mix * sigma_mix )
	923
	924	IF ( qvapinmix .GT. qsat(i) ) THEN
	925	!--If the mixed air is supersaturated, we condense the subsaturated
	926	!--region which was mixed.
	927	dcf_mix_sub = subfra_mix
	928	dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
	929	dqvc_mix_sub = dcf_mix_sub * qsat(i)
	930	ELSE
	931	!--Else, we sublimate the cloud which was mixed.
	932	dcf_mix_sub = - sigma_mix * cldfra_mix
	933	dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
	934	dqvc_mix_sub = dcf_mix_sub * qsat(i)
	935	ENDIF
	936	ENDIF ! ok_unadjusted_clouds
	937	ENDIF ! subfra .GT. eps
	938
	939	!--We then mix the supersaturated sky (issrfra_mix) and the cloud,
	940	!--for which the mixed air is always supersatured, therefore
	941	!--the cloud necessarily expands
	942	IF ( issrfra(i) .GT. eps ) THEN
	943
	944	IF ( ok_unadjusted_clouds ) THEN
	945	!--The ice supersaturated air is simply added to the
	946	!--cloud, and supersaturated vapor will be deposited on the
	947	!--cloud ice crystals by the deposition process in the
	948	!--following timestep
	949	dcf_mix_issr = issrfra_mix
	950	dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
	951	dqvc_mix_issr = dqt_mix_issr
	952	ELSE
	953	!--In this case, the additionnal vapor condenses
	954	dcf_mix_issr = issrfra_mix
	955	dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
	956	dqvc_mix_issr = dcf_mix_issr * qsat(i)
	957	ENDIF ! ok_unadjusted_clouds
	958
	959
	960	ENDIF ! issrfra .GT. eps
	961
	962	!--Sum up the tendencies from subsaturated sky and supersaturated sky
	963	dcf_mix(i) = dcf_mix_sub + dcf_mix_issr
	964	dqt_mix = dqt_mix_sub + dqt_mix_issr
	965	dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
	966	dqi_mix(i) = dqt_mix - dqvc_mix(i)
	967
	968	!--Add tendencies
	969	issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
	970	qissr(i) = MAX(0., qissr(i) - dqt_mix_issr)
	971	cldfra(i) = MAX(0., MIN(1. - dcf_con(i), cldfra(i) + dcf_mix(i)))
	972	qcld(i) = MAX(0., MIN(qtot(i) - dqi_con(i) - dqvc_con(i), qcld(i) + dqt_mix))
	973	qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
	974
	975	ENDIF ! ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
	976
	977	!--Finally, we add the tendencies of condensation
	978	cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))
	979	qcld(i) = MIN(qtot(i), qcld(i) + dqvc_con(i) + dqi_con(i))
	980	qvc(i) = MIN(qcld(i), qvc(i) + dqvc_con(i))
	981
	982
	983	!----------------------------------------
	984	!-- CONTRAILS AND AVIATION --
	985	!----------------------------------------
	986
	987	!--Add a source of cirrus from aviation contrails
	988	!IF ( ok_plane_contrail ) THEN
	989	! dcf_avi(i) = 0.
	990	! dqi_avi(i) = 0.
	991	! dqvc_avi(i) = 0.
	992	! ! TODO implement ok_unadjusted_clouds
	993	! IF ( issrfra(i) .GT. eps ) THEN
	994	! contrail_fra = MIN(1., flight_m(i,k) * dtime * contrail_cross_section / V_cell)
	995	! dcf_avi(i) = issrfra(i) * contrail_fra
	996	! dqt_avi = dcf_avi(i) * qissr(i) / issrfra(i)
	997	! dqvc_avi(i) = qsat(i) * dcf_avi(i)
	998	!
	999	! !--Add tendencies
	1000	! cldfra(i) = cldfra(i) + dcf_avi(i)
	1001	! issrfra(i) = issrfra(i) - dcf_avi(i)
	1002	! qcld(i) = qcld(i) + dqt_avi
	1003	! qvc(i) = qvc(i) + dqvc_avi(i)
	1004	! qissr(i) = qissr(i) - dqt_avi
	1005
	1006	! !--Diagnostics
	1007	! dqi_avi(i) = dqt_avi - qsat(i) * dcf_avi(i)
	1008	! ENDIF
	1009	! dcf_avi(i) = dcf_avi(i) / dtime
	1010	! dqi_avi(i) = dqi_avi(i) / dtime
	1011	! dqvc_avi(i) = dqvc_avi(i) / dtime
	1012	!ENDIF
	1013
	1014
	1015
	1016	!-------------------------------------------
	1017	!-- FINAL BARRIERS AND OUTPUTS --
	1018	!-------------------------------------------
	1019
	1020	IF ( cldfra(i) .LT. eps ) THEN
	1021	!--If the cloud is too small, it is sublimated.
	1022	cldfra(i) = 0.
	1023	qcld(i) = 0.
	1024	qvc(i) = 0.
	1025	qincld(i) = qsat(i)
	1026	ELSE
	1027	qincld(i) = qcld(i) / cldfra(i)
	1028	ENDIF ! cldfra .LT. eps
	1029
	1030	!--Diagnostics
	1031	dcf_sub(i) = dcf_sub(i) / dtime
	1032	dcf_con(i) = dcf_con(i) / dtime
	1033	dcf_mix(i) = dcf_mix(i) / dtime
	1034	dqi_adj(i) = dqi_adj(i) / dtime
	1035	dqi_sub(i) = dqi_sub(i) / dtime
	1036	dqi_con(i) = dqi_con(i) / dtime
	1037	dqi_mix(i) = dqi_mix(i) / dtime
	1038	dqvc_adj(i) = dqvc_adj(i) / dtime
	1039	dqvc_sub(i) = dqvc_sub(i) / dtime
	1040	dqvc_con(i) = dqvc_con(i) / dtime
	1041	dqvc_mix(i) = dqvc_mix(i) / dtime
	1042
	1043	ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
	1044
	1045	ENDIF ! end keepgoing
	1046
	1047	ENDDO ! end loop on i
	1048
	1049	END SUBROUTINE condensation_ice_supersat
	1050	!**********************************************************************************
	1051
	1052	END MODULE lmdz_lscp_condensation

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

Download in other formats:

Original Format