Context Navigation

lscp_mod.F90 @ 4190

Last change on this file since 4190 was 4013, checked in by lguez, 3 years ago
Sync latest trunk changes to Ocean_skin
File size: 51.5 KB

Rev	Line
[3999]	1	MODULE LSCP_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
	7	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	8	SUBROUTINE LSCP(dtime,paprs,pplay,t,q,ptconv,ratqs, &
	9	d_t, d_q, d_ql, d_qi, rneb, radliq, radicefrac, &
	10	rain, snow, &
	11	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
	12	frac_impa, frac_nucl, beta, &
	13	prfl, psfl, rhcl, zqta, fraca, &
	14	ztv, zpspsk, ztla, zthl, iflag_cld_th, &
	15	iflag_ice_thermo)
	16
	17	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	18	! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
	19	! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
	20	!--------------------------------------------------------------------------------
	21	! Date: 01/2021
	22	!--------------------------------------------------------------------------------
	23	! Aim: Large Scale Clouds and Precipitation (LSCP)
	24	!
	25	! This code is a new version of the fisrtilp.F90 routine, which is itself a
	26	! fusion of 'first' (superrsaturation physics, P. LeVan K. Laval)
	27	! and 'ilp' (il pleut, L. Li)
	28	!
	29	! Compared to fisrtilp, LSCP
	30	! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
	31	! -> consider always precipitation thermalisation (fl_cor_ebil>0)
	32	! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
	33	! -> option "oldbug" by JYG has been removed
	34	! -> iflag_t_glace >0 always
	35	! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
	36	! -> rectangular distribution from L. Li no longer available
	37	! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
	38	!--------------------------------------------------------------------------------
	39	! References:
	40	!
	41	! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
	42	! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
	43	! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
	44	! -------------------------------------------------------------------------------
	45	! Code structure:
	46	!
	47	! P0> Thermalization of the precipitation coming from the overlying layer
	48	! P1> Evaporation of the precipitation (falling from the k+1 level)
	49	! P2> Cloud formation (at the k level)
	50	! P2.A.0> Cloud properties calculation from a rectangular pdf
	51	! P2.A.1> With the new PDFs, calculation of cloud properties using the inital
	52	! values of T and Q
	53	! P2.A.2> Coupling between condensed water and temperature
	54	! P2.A.3> Calculation of final quantities associated with cloud formation
	55	! P2.B> 'All or nothing' cloud
	56	! P2.C> Release of Latent heat after cloud formation
	57	! P3> Autoconversion to precipitation (k-level)
	58	! P4> Wet scavenging
	59	!------------------------------------------------------------------------------
	60	! Some preliminary comments (JBM) :
	61	!
	62	! The cloud water that the radiation scheme sees is not the same that the cloud
	63	! water used in the physics and the dynamics
	64	!
	65	! During the autoconversion to precipitation (P3 step), radliq (cloud water used
	66	! by the radiation scheme) is calculated as an average of the water that remains
	67	! in the cloud during the precipitation and not the water remaining at the end
	68	! of the time step. The latter is used in the rest of the physics and advected
	69	! by the dynamics.
	70	!
	71	! In summary:
	72	!
	73	! Radiation:
	74	! xflwc(newmicro)+xfiwc(newmicro) =
	75	! cldliq(physiq)=radliq(fisrt)=lwcon(nc)+iwcon(nc)
	76	!
	77	! Physics/Dynamics:
	78	! ql_seri(physiq)+qs_seri(physiq)=ocond(nc)
	79	!
	80	! Notetheless, be aware of:
	81	!
	82	! radliq(fisrt) .NE. ocond(nc)
	83	! i.e.:
	84	! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
	85	! (which is not trivial)
	86	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	87
	88	USE dimphy
	89	USE print_control_mod, ONLY: prt_level, lunout
	90	USE cloudth_mod
	91	USE ioipsl_getin_p_mod, ONLY : getin_p
	92	USE phys_local_var_mod, ONLY: ql_seri,qs_seri
	93	USE phys_local_var_mod, ONLY: rneblsvol
	94	USE lscp_tools_mod, ONLY : CALC_QSAT_ECMWF, ICEFRAC_LSCP, CALC_GAMMASAT, FALLICE_VELOCITY
	95
	96
	97	IMPLICIT NONE
	98
	99
	100
	101	!===============================================================================
	102	! VARIABLE DECLARATION
	103	!===============================================================================
	104
	105	include "YOMCST.h"
	106	include "YOETHF.h"
	107	include "FCTTRE.h"
	108	include "fisrtilp.h"
	109	include "nuage.h"
	110
	111
	112	! INPUT VARIABLES:
	113	!-----------------
	114
	115	REAL, INTENT(IN) :: dtime ! time step [s]
	116	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa]
	117	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa]
	118	REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
	119	REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! specific humidity [kg/kg]
	120	INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds
	121	INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
	122	! CR: if iflag_ice_thermo=2, only convection is active
	123	LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
	124
	125
	126	!Inputs associated with thermal plumes
	127
	128	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! virtual potential temperature [K]
	129	REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg]
	130	REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-]
	131	REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk ! exner potential (p/100000)**(R/cp)
	132	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! liquid temperature within thermals [K]
	133	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! liquid potential temperature [K]
	134
	135	! INPUT/OUTPUT variables
	136	!------------------------
	137
	138	REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! function of pressure that sets the large-scale
	139	! cloud PDF (sigma=ratqs*qt)
	140
	141	! OUTPUT variables
	142	!-----------------
	143
	144	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K]
	145	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg]
	146	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg]
	147	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg]
	148	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-]
	149	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radliq ! condensed water used in the radiation scheme [kg/kg]
	150	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme
	151	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-]
	152	REAL, DIMENSION(klon), INTENT(OUT) :: rain ! large-scale rainfall [kg/s/m2]
	153	REAL, DIMENSION(klon), INTENT(OUT) :: snow ! large-scale snowfall [kg/s/m2]
	154	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2]
	155	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2]
	156
	157
	158	! cofficients of scavenging fraction (for off-line)
	159
	160	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_nucl ! scavenging fraction due tu nucleation [-]
	161	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_1nucl ! scavenging fraction due tu nucleation with a -1 factor [-]
	162	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_impa ! scavening fraction due to impaction [-]
	163
	164	! fraction of aerosol scavenging through impaction and nucleation (for on-line)
	165
	166	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-]
	167	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-]
	168
	169
	170
	171
	172	! PROGRAM PARAMETERS:
	173	!--------------------
	174
	175
	176	REAL, SAVE :: seuil_neb=0.001 ! cloud fraction threshold: a cloud really exists when exceeded
	177	!$OMP THREADPRIVATE(seuil_neb)
	178
	179	INTEGER, SAVE :: ninter=5 ! number of iterations to calculate autoconversion to precipitation
	180	!$OMP THREADPRIVATE(ninter)
	181
	182	INTEGER,SAVE :: iflag_evap_prec=1 ! precipitation evap. flag. 0: nothing, 1: "old way", 2: Max cloud fraction above to calculate the max of reevaporation
	183	! 4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only
	184	!$OMP THREADPRIVATE(iflag_evap_prec)
	185
	186
	187	REAL t_coup ! temperature threshold which determines the phase
	188	PARAMETER (t_coup=234.0) ! for which the saturation vapor pressure is calculated
	189
	190	REAL DDT0 ! iteration parameter
	191	PARAMETER (DDT0=.01)
	192
	193	REAL ztfondue ! parameter to calculate melting fraction of precipitation
	194	PARAMETER (ztfondue=278.15)
	195
	196	REAL, SAVE :: rain_int_min=0.001 ! Minimum local rain intensity [mm/s] before the decrease in associates precipitation fraction
	197	!$OMP THREADPRIVATE(rain_int_min)
	198
	199
	200
	201	! LOCAL VARIABLES:
	202	!----------------
	203
	204
	205	REAL qsl, qsi
	206	REAL zct, zcl
	207	REAL ctot(klon,klev)
	208	REAL ctot_vol(klon,klev)
	209	INTEGER mpc_bl_points(klon,klev)
	210	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
	211	REAL zdqsdT_raw(klon)
	212	REAL gammasat(klon),dgammasatdt(klon) ! coefficient to make cold condensation at the correct RH and derivative wrt T
	213	REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
	214	REAL cste
	215	REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
	216	REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
	217	REAL erf
	218	REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon)
	219	REAL zrfl(klon), zrfln(klon), zqev, zqevt
	220	REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
	221	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
	222	REAL zoliqp(klon), zoliqi(klon)
	223	REAL zt(klon)
	224	REAL zdz(klon),zrho(klon),ztot, zrhol(klon)
	225	REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
	226	REAL zmelt,zpluie,zice
	227	REAL dzfice(klon)
	228	REAL zsolid,qtot,dqsl,dqsi
	229	REAL smallestreal
	230	! Variables for Bergeron process
	231	REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
	232	REAL zqpreci(klon), zqprecl(klon)
	233	! Variables precipitation energy conservation
	234	REAL zmqc(klon)
	235	! Variables for tracers
	236	! temporary: alpha parameter for scavenging
	237	! 4 possible scavenging processes
	238	REAL a_tr_sca(4)
	239	save a_tr_sca
	240	!$OMP THREADPRIVATE(a_tr_sca)
	241	REAL zalpha_tr
	242	REAL zfrac_lessi
	243	REAL zprec_cond(klon)
	244	REAL beta(klon,klev) ! conversion rate of condensed water
	245	REAL zmair(klon), zcpair, zcpeau
	246	REAL zlh_solid(klon), zm_solid ! for liquid -> solid conversion
	247	REAL zrflclr(klon), zrflcld(klon)
	248	REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)
	249	REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon)
	250	REAL ziflclr(klon), ziflcld(klon)
	251	REAL znebprecipclr(klon), znebprecipcld(klon)
	252	REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
	253	REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
	254	REAL velo(klon)
	255	REAL qlmpc, qimpc, rnebmpc
	256	REAL radliqi(klon,klev)
	257
	258	INTEGER i, k, n, kk
	259	INTEGER n_i(klon), iter
	260	INTEGER ncoreczq
	261	INTEGER,SAVE :: itap=0
	262	!$OMP THREADPRIVATE(itap)
	263
	264
	265	LOGICAL lognormale(klon)
	266	LOGICAL convergence(klon)
	267	LOGICAL appel1er
	268	SAVE appel1er
	269	!$OMP THREADPRIVATE(appel1er)
	270
	271
	272	CHARACTER (len = 20) :: modname = 'lscp'
	273	CHARACTER (len = 80) :: abort_message
	274
	275
	276	DATA appel1er /.TRUE./
	277
	278
	279
	280
	281	!===============================================================================
	282	! INITIALISATION
	283	!===============================================================================
	284
	285	! Few initial checksS
	286
	287	IF (iflag_t_glace.EQ.0) THEN
	288	abort_message = 'lscp cannot be used if iflag_t_glace=0'
	289	CALL abort_physic(modname,abort_message,1)
	290	ENDIF
	291
	292	IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
	293	abort_message = 'lscp cannot be used without ice thermodynamics'
	294	CALL abort_physic(modname,abort_message,1)
	295	ENDIF
	296
	297	IF (.NOT.thermcep) THEN
	298	abort_message = 'lscp cannot be used when thermcep=false'
	299	CALL abort_physic(modname,abort_message,1)
	300	ENDIF
	301
	302	IF (iflag_fisrtilp_qsat .LT. 0) THEN
	303	abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
	304	CALL abort_physic(modname,abort_message,1)
	305	ENDIF
	306
	307
	308	! Few initialisations
	309
	310	itap=itap+1
	311	znebprecip(:)=0.
	312	ctot_vol(1:klon,1:klev)=0.0
	313	rneblsvol(1:klon,1:klev)=0.0
	314	smallestreal=1.e-9
	315	znebprecipclr(:)=0.
	316	znebprecipcld(:)=0.
	317	mpc_bl_points(:,:)=0
	318
	319	IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
	320
	321	IF (appel1er) THEN
	322	CALL getin_p('ninter',ninter)
	323	CALL getin_p('iflag_evap_prec',iflag_evap_prec)
	324	CALL getin_p('seuil_neb',seuil_neb)
	325	CALL getin_p('rain_int_min',rain_int_min)
	326	WRITE(lunout,*) 'lscp, ninter:', ninter
	327	WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec
	328	WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb
	329	WRITE(lunout,*) 'lscp, rain_int_min:', rain_int_min
	330
	331	! check for precipitation sub-time steps
	332	IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
	333	WRITE(lunout,*) 'lscp: it is not expected, see Z.X.Li', dtime
	334	WRITE(lunout,*) 'I would prefer a 6 min sub-timestep'
	335	ENDIF
	336
	337	!AA Temporary initialisation
	338	a_tr_sca(1) = -0.5
	339	a_tr_sca(2) = -0.5
	340	a_tr_sca(3) = -0.5
	341	a_tr_sca(4) = -0.5
	342
	343	!AA Set scavenged fractions to 1
	344	DO k = 1, klev
	345	DO i = 1, klon
	346	pfrac_nucl(i,k)=1.
	347	pfrac_1nucl(i,k)=1.
	348	pfrac_impa(i,k)=1.
	349	beta(i,k)=0.
	350	ENDDO
	351	ENDDO
	352
	353	appel1er = .FALSE.
	354
	355	ENDIF !(appel1er)
	356
	357	! AA for 'safety' reasons
	358	zalpha_tr = 0.
	359	zfrac_lessi = 0.
	360
	361
	362	! Initialisation of output variables (JYG):
	363	prfl(:,:) = 0.0
	364	psfl(:,:) = 0.0
	365
	366	d_t(:,:) = 0.0
	367	d_q(:,:) = 0.0
	368	d_ql(:,:) = 0.0
	369	d_qi(:,:) = 0.0
	370	rneb(:,:) = 0.0
	371	radliq(:,:) = 0.0
	372	radicefrac(:,:) = 0.0
	373	frac_nucl(:,:) = 1.
	374	frac_impa(:,:) = 1.
	375
	376	rain(:) = 0.0
	377	snow(:) = 0.0
	378	zoliq(:)=0.
	379	zrfl(:) = 0.0
	380	zifl(:) = 0.0
	381	zneb(:) = seuil_neb
	382	zrflclr(:) = 0.0
	383	ziflclr(:) = 0.0
	384	zrflcld(:) = 0.0
	385	ziflcld(:) = 0.0
	386	tot_zneb(:) = 0.0
	387	tot_znebn(:) = 0.0
	388	d_tot_zneb(:) = 0.0
	389
	390
	391	!===============================================================================
	392	! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
	393	!===============================================================================
	394
	395	ncoreczq=0
	396
	397	DO k = klev, 1, -1
	398
	399
	400	! Initialisation temperature and specific humidity
	401	DO i = 1, klon
	402	zt(i)=t(i,k)
	403	zq(i)=q(i,k)
	404	ENDDO
	405
	406	! --------------------------------------------------------------------
	407	! P0> Thermalization of precipitation falling from the overlying layer
	408	! --------------------------------------------------------------------
	409	! Computes air temperature variation due to latent heat transported by
	410	! precipitation. Precipitation is then thermalized with the air in the
	411	! layer.
	412	! The precipitation should remain thermalized throughout the different
	413	! thermodynamical transformations. The corresponding water mass should
	414	! be added when calculating the layer's enthalpy change with
	415	! temperature
	416	! ---------------------------------------------------------------------
	417
	418	IF (k.LE.klevm1) THEN
	419
	420	DO i = 1, klon
	421
	422	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
	423	! no condensed water so cp=cp(vapor+dry air)
	424	! RVTMP2=rcpv/rcpd-1
	425	zcpair=RCPD(1.0+RVTMP2zq(i))
	426	zcpeau=RCPD*RVTMP2
	427
	428	! zmqc: precipitation mass that has to be thermalized with
	429	! layer's air so that precipitation at the ground has the
	430	! same temperature as the lowermost layer
	431	zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
	432	! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
	433	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zmqc(i)zcpeau + zcpair*zt(i) ) &
	434	/ (zcpair + zmqc(i)*zcpeau)
	435
	436	ENDDO
	437
	438	ELSE
	439
	440	DO i = 1, klon
	441	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
	442	zmqc(i) = 0.
	443	ENDDO
	444
	445	ENDIF
	446
	447	! --------------------------------------------------------------------
	448	! P1> Precipitation evaporation/sublimation
	449	! --------------------------------------------------------------------
	450	! A part of the precipitation coming from above is evaporated/sublimated.
	451	!
	452	! --------------------------------------------------------------------
	453
	454
	455	IF (iflag_evap_prec.GE.1) THEN
	456
	457
	458	DO i = 1, klon
	459
	460	! if precipitation
	461	IF (zrfl(i)+zifl(i).GT.0.) THEN
	462
	463	CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
	464
	465	! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
	466	IF (iflag_evap_prec.EQ.4) THEN
	467	zrfl(i) = zrflclr(i)
	468	zifl(i) = ziflclr(i)
	469	ENDIF
	470
	471	IF (iflag_evap_prec.EQ.1) THEN
	472	znebprecip(i)=zneb(i)
	473	ELSE
	474	znebprecip(i)=MAX(zneb(i),znebprecip(i))
	475	ENDIF
	476
	477
	478	IF (iflag_evap_prec.EQ.4) THEN
	479	! Max evaporation not to saturate the whole mesh
	480	zqev0 = MAX(0.0, zqs(i)-zq(i))
	481	ELSE
	482	! Max evap not to saturate the fraction below the cloud
	483	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
	484	ENDIF
	485
	486	! A. JAM
	487	! We consider separately qsatice and qstal
	488	! qsat wrt liquid phase according to thermcep
	489	CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,1,.false.,qsl,dqsl)
	490
	491
	492	! Evaporation of liquid precipitation coming from above
	493	! dP/dz=beta(1-q/qsat)sqrt(P)
	494	! formula from Sundquist 1988, Klemp & Wilhemson 1978
	495	! LTP: evaporation only in the clear sky part (iflag_evap_prec=4)
	496
	497	IF (iflag_evap_prec.EQ.3) THEN
	498	zqevt = znebprecip(i)coef_eva(1.0-zq(i)/qsl) &
	499	*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
	500	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	501	ELSE IF (iflag_evap_prec.EQ.4) THEN
	502	zqevt = znebprecipclr(i)coef_eva(1.0-zq(i)/qsl) &
	503	*SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
	504	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	505	ELSE
	506	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
	507	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	508	ENDIF
	509
	510
	511	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
	512	RGdtime/(paprs(i,k)-paprs(i,k+1))
	513
	514	! qsat wrt ice according to thermcep
	515	CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,2,.false.,qsi,dqsi)
	516
	517	! sublimation of the solid precipitation coming from above
	518	IF (iflag_evap_prec.EQ.3) THEN
	519	zqevti = znebprecip(i)coef_eva(1.0-zq(i)/qsi) &
	520	*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
	521	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	522	ELSE IF (iflag_evap_prec.EQ.4) THEN
	523	zqevti = znebprecipclr(i)coef_eva(1.0-zq(i)/qsi) &
	524	*SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
	525	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	526	ELSE
	527	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
	528	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	529	ENDIF
	530
	531	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
	532	RGdtime/(paprs(i,k)-paprs(i,k+1))
	533
	534	! A. JAM
	535	! Evaporation limit: we ensure that the layer's fraction below
	536	! the cloud or the whole mesh (depending on iflag_evap_prec)
	537	! does not reach saturation. In this case, we
	538	! redistribute zqev0 conserving the ratio liquid/ice
	539
	540	IF (zqevt+zqevti.GT.zqev0) THEN
	541	zqev=zqev0*zqevt/(zqevt+zqevti)
	542	zqevi=zqev0*zqevti/(zqevt+zqevti)
	543	ELSE
	544	zqev=zqevt
	545	zqevi=zqevti
	546	ENDIF
	547
	548
	549	! EV: agrees with JL's comments below so correct and comment
	550	! JLD: I do not understand the lines below. We distribute the liquid
	551	! and solid parts of the precipitation even if the layer does not get
	552	! saturated. To my opinion, we should all replace with:
	553	! zqev=zqevt
	554	! zqevi=zqevti
	555	! IF (zqevt+zqevti.GT.0.) THEN
	556	! zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
	557	! zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
	558	! ELSE
	559	! zqev=0.
	560	! zqevi=0.
	561	! ENDIF
	562	! ENDIF
	563
	564	! New solid and liquid precipitation fluxes
	565	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
	566	/RG/dtime)
	567	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
	568	/RG/dtime)
	569
	570	! vapor, temperature, precip fluxes update
	571	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	572	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
	573	! precip thermalization
	574	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	575	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
	576	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	577	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	578	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
	579	+ (zifln(i)-zifl(i)) &
	580	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	581	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
	582
	583	! New values of liquid and solid precipitation
	584	zrfl(i) = zrfln(i)
	585	zifl(i) = zifln(i)
	586
	587	IF (iflag_evap_prec.EQ.4) THEN
	588	zrflclr(i) = zrfl(i)
	589	ziflclr(i) = zifl(i)
	590	IF(zrflclr(i) + ziflclr(i).LE.0) THEN
	591	znebprecipclr(i) = 0.
	592	ENDIF
	593	zrfl(i) = zrflclr(i) + zrflcld(i)
	594	zifl(i) = ziflclr(i) + ziflcld(i)
	595	ENDIF
	596
	597
	598	! CR Be careful: ice melting has been moved
	599	zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! JYG
	600	! precip fraction that is melted
	601	zmelt = MIN(MAX(zmelt,0.),1.)
	602	! Ice melting
	603	IF (iflag_evap_prec.EQ.4) THEN
	604	zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
	605	zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
	606	zrfl(i)=zrflclr(i)+zrflcld(i)
	607	ELSE
	608	zrfl(i)=zrfl(i)+zmelt*zifl(i)
	609	ENDIF
	610
	611	! Latent heat of melting with precipitation thermalization
	612	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
	613	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
	614
	615	IF (iflag_evap_prec.EQ.4) THEN
	616	ziflclr(i)=ziflclr(i)*(1.-zmelt)
	617	ziflcld(i)=ziflcld(i)*(1.-zmelt)
	618	zifl(i)=ziflclr(i)+ziflcld(i)
	619	ELSE
	620	zifl(i)=zifl(i)*(1.-zmelt)
	621	ENDIF
	622
	623	ELSE
	624	! if no precip, we reinitialize the cloud fraction used for the precip to 0
	625	znebprecip(i)=0.
	626
	627	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	628
	629	ENDDO
	630
	631	ENDIF ! (iflag_evap_prec>=1)
	632
	633	! --------------------------------------------------------------------
	634	! End precip evaporation
	635	! --------------------------------------------------------------------
	636
	637
	638
	639
	640	! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
	641	!-------------------------------------------------------
	642
	643	DO i = 1, klon
	644	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
	645	qtot=zq(i)+zmqc(i)
	646	CALL CALC_QSAT_ECMWF(zt(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
	647	zdqsdT_raw(i) = zdqs(i)* &
	648	& RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTTzdelta)
	649
	650	IF (zq(i) .LT. 1.e-15) THEN
	651	ncoreczq=ncoreczq+1
	652	zq(i)=1.e-15
	653	ENDIF
	654
	655	ENDDO
	656
	657
	658	! --------------------------------------------------------------------
	659	! P2> Cloud formation
	660	!---------------------------------------------------------------------
	661	!
	662	! Unlike fisrtilp, we always assume a 'fractional cloud' approach
	663	! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
	664	! is prescribed and depends on large scale variables and boundary layer
	665	! properties)
	666	! The decrease in condensed part due tu latent heating is taken into
	667	! account
	668	! -------------------------------------------------------------------
	669
	670
	671	! P2.1> With the PDFs (log-normal, bigaussian)
	672	! cloud propertues calculation with the initial values of t and q
	673	! ----------------------------------------------------------------
	674
	675	! initialise gammasat and qincloud_mpc
	676	gammasat(:)=1.
	677	qincloud_mpc(:)=0.
	678
	679
	680
	681	IF (iflag_cld_th.GE.5) THEN
	682
	683	IF (iflag_mpc_bl .LT. 1) THEN
	684
	685	IF (iflag_cloudth_vert.LE.2) THEN
	686
	687	CALL cloudth(klon,klev,k,ztv, &
	688	zq,zqta,fraca, &
	689	qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
	690	ratqs,zqs,t)
	691
	692	ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
	693
	694	CALL cloudth_v3(klon,klev,k,ztv, &
	695	zq,zqta,fraca, &
	696	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
	697	ratqs,zqs,t)
	698
	699
	700	!Jean Jouhaud's version, Decembre 2018
	701	!-------------------------------------
	702
	703	ELSEIF (iflag_cloudth_vert.EQ.6) THEN
	704
	705	CALL cloudth_v6(klon,klev,k,ztv, &
	706	zq,zqta,fraca, &
	707	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
	708	ratqs,zqs,t)
	709
	710	ENDIF
	711
	712	ELSE
	713	! cloudth_v3 + cold microphysical considerations
	714	! + stationary mixed-phase cloud model
	715
	716	CALL cloudth_mpc(klon,klev,k,mpc_bl_points, &
	717	t,ztv,zq,zqta,fraca, zpspsk, &
	718	paprs,pplay,ztla,zthl,ratqs,zqs,psfl, &
	719	qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol)
	720
	721	ENDIF ! iflag_mpc_bl
	722
	723	DO i=1,klon
	724	rneb(i,k)=ctot(i,k)
	725	rneblsvol(i,k)=ctot_vol(i,k)
	726	zqn(i)=qcloud(i)
	727	ENDDO
	728
	729	ENDIF
	730
	731	IF (iflag_cld_th .LE. 4) THEN
	732
	733	! lognormal
	734	lognormale = .TRUE.
	735
	736	ELSEIF (iflag_cld_th .GE. 6) THEN
	737
	738	! lognormal distribution when no thermals
	739	lognormale = fraca(:,k) < 1e-10
	740
	741	ELSE
	742	! When iflag_cld_th=5, we always assume
	743	! bi-gaussian distribution
	744	lognormale = .FALSE.
	745
	746	ENDIF
	747
	748
	749	DT(:) = 0.
	750	n_i(:)=0
	751	Tbef(:)=zt(:)
	752	qlbef(:)=0.
	753
	754	! Treatment of non-boundary layer clouds (lognormale)
	755	! condensation with qsat(T) variation (adaptation)
	756	! Iterative Loop to converge towards qsat
	757
	758
	759	DO iter=1,iflag_fisrtilp_qsat+1
	760
	761	DO i=1,klon
	762
	763	! convergence = .true. since convergence is not satisfied
	764	convergence(i)=ABS(DT(i)).GT.DDT0
	765
	766	IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
	767
	768	! if not convergence:
	769
	770	! P2.2.1> cloud fraction and condensed water mass calculation
	771	! Calculated variables:
	772	! rneb : cloud fraction
	773	! zqn : total water within the cloud
	774	! zcond: mean condensed water within the mesh
	775	! rhcl: clear-sky relative humidity
	776	!---------------------------------------------------------------
	777
	778	! new temperature:
	779	Tbef(i)=Tbef(i)+DT(i)
	780
	781
	782	! Rneb, qzn and zcond for lognormal PDFs
	783	qtot=zq(i)+zmqc(i)
	784	CALL CALC_QSAT_ECMWF(Tbef(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
	785	CALL CALC_GAMMASAT(Tbef(i),qtot,pplay(i,k),gammasat(i),dgammasatdt(i))
	786
	787	! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
	788	zdqs(i) = gammasat(i)zdqs(i)+zqs(i)dgammasatdt(i)
	789
	790	zpdf_sig(i)=ratqs(i,k)*zq(i)
	791	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
	792	zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
	793	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
	794	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
	795	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
	796	zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
	797	zpdf_e1(i)=1.-erf(zpdf_e1(i))
	798	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
	799	zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
	800	zpdf_e2(i)=1.-erf(zpdf_e2(i))
	801
	802	IF (zpdf_e1(i).LT.1.e-10) THEN
	803	rneb(i,k)=0.
	804	zqn(i)=gammasat(i)*zqs(i)
	805	ELSE
	806	rneb(i,k)=0.5*zpdf_e1(i)
	807	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
	808	ENDIF
	809
	810	! If vertical heterogeneity, change fraction by volume as well
	811	IF (iflag_cloudth_vert.GE.3) THEN
	812	ctot_vol(i,k)=rneb(i,k)
	813	rneblsvol(i,k)=ctot_vol(i,k)
	814	ENDIF
	815
	816
	817	! P2.2.2> Approximative calculation of temperature variation DT
	818	! due to condensation.
	819	! Calculated variables:
	820	! dT : temperature change due to condensation
	821	! ---------------------------------------------------------------
	822
	823	! EV: calculation of icefrac in one sole function
	824	CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
	825
	826
	827	IF (zfice(i).LT.1) THEN
	828	cste=RLVTT
	829	ELSE
	830	cste=RLSTT
	831	ENDIF
	832
	833	qlbef(i)=max(0.,zqn(i)-zqs(i))
	834	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
	835	& +zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
	836	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
	837	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
	838	& qlbef(i)dzfice(i)
	839	DT(i)=num/denom
	840	n_i(i)=n_i(i)+1
	841
	842	ENDIF ! end convergence
	843
	844	ENDDO ! end loop on i
	845
	846	ENDDO ! iter=1,iflag_fisrtilp_qsat+1
	847
	848
	849
	850	! P2.3> Final quantities calculation
	851	! Calculated variables:
	852	! rneb : cloud fraction
	853	! zcond: mean condensed water in the mesh
	854	! zqn : mean water vapor in the mesh
	855	! zt : temperature
	856	! rhcl : clear-sky relative humidity
	857	! ----------------------------------------------------------------
	858
	859
	860
	861	! Water vapor update, Phase determination and subsequent latent heat exchange
	862
	863	! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
	864	CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
	865
	866	DO i=1, klon
	867
	868
	869	IF (mpc_bl_points(i,k) .GT. 0) THEN
	870	zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
	871	! following line is very strange and probably wrong
	872	rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
	873	! water vapor update and partition function if thermals
	874	zq(i) = zq(i) - zcond(i)
	875	zfice(i)=icefrac_mpc(i,k)
	876
	877	ELSE
	878
	879	! Checks on rneb, rhcl and zqn
	880	IF (rneb(i,k) .LE. 0.0) THEN
	881	zqn(i) = 0.0
	882	rneb(i,k) = 0.0
	883	zcond(i) = 0.0
	884	rhcl(i,k)=zq(i)/zqs(i)
	885	ELSE IF (rneb(i,k) .GE. 1.0) THEN
	886	zqn(i) = zq(i)
	887	rneb(i,k) = 1.0
	888	zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))
	889	rhcl(i,k)=1.0
	890	ELSE
	891	zcond(i) = MAX(0.0,zqn(i)-gammasat(i)zqs(i))rneb(i,k)
	892	! following line is very strange and probably wrong:
	893	rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
	894	ENDIF
	895
	896
	897	! water vapor update
	898	zq(i) = zq(i) - zcond(i)
	899
	900	ENDIF
	901
	902	! temperature update due to phase change
	903	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
	904	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
	905	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
	906	ENDDO
	907
	908
	909	! If vertical heterogeneity, change volume fraction
	910	IF (iflag_cloudth_vert .GE. 3) THEN
	911	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
	912	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
	913	ENDIF
	914
	915	! ----------------------------------------------------------------
	916	! P3> Precipitation formation and calculation of condensed
	917	! water seen by the radiation scheme
	918	! ----------------------------------------------------------------
	919
	920	! LTP:
	921	IF (iflag_evap_prec==4) THEN
	922
	923	!Partitionning between precipitation coming from clouds and that coming from CS
	924
	925	!0) Calculate tot_zneb, total cloud fraction above the cloud
	926	!assuming a maximum-random overlap (voir Jakob and Klein, 2000)
	927
	928	DO i=1, klon
	929	tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
	930	/(1-min(zneb(i),1-smallestreal))
	931	d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
	932	tot_zneb(i) = tot_znebn(i)
	933
	934
	935	!1) Cloudy to clear air
	936	d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
	937	IF (znebprecipcld(i) .GT. 0) THEN
	938	d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
	939	d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
	940	ELSE
	941	d_zrfl_cld_clr(i) = 0.
	942	d_zifl_cld_clr(i) = 0.
	943	ENDIF
	944
	945	!2) Clear to cloudy air
	946	d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) &
	947	- d_tot_zneb(i) - zneb(i)))
	948	IF (znebprecipclr(i) .GT. 0) THEN
	949	d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
	950	d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
	951	ELSE
	952	d_zrfl_clr_cld(i) = 0.
	953	d_zifl_clr_cld(i) = 0.
	954	ENDIF
	955
	956	!Update variables
	957	znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
	958	znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
	959	zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
	960	ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
	961	zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
	962	ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
	963
	964	ENDDO
	965	ENDIF
	966
	967	! Initialisation of zoliq and radliq variables
	968
	969	DO i = 1, klon
	970	IF (rneb(i,k).GT.0.0) THEN
	971	zoliq(i) = zcond(i)
	972	zrho(i) = pplay(i,k) / zt(i) / RD
	973	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
	974
	975	zneb(i) = MAX(rneb(i,k), seuil_neb)
	976	radliq(i,k) = zoliq(i)/REAL(ninter+1)
	977	radicefrac(i,k) = zfice(i)
	978	radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
	979	ENDIF
	980	ENDDO
	981
	982	! -------------------------------------------------------------------------------
	983	! Iterations to calculate a "mean" cloud water content during the precipitation
	984	! Comment: it is not the remaining water after precipitation but a mean
	985	! remaining water in the cloud during the time step that is seen by the radiation
	986	! -------------------------------------------------------------------------------
	987
	988
	989	DO n = 1, ninter
	990
	991	DO i=1,klon
	992	IF (rneb(i,k).GT.0.0) THEN
	993	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	994	ENDIF
	995	ENDDO
	996
	997	CALL FALLICE_VELOCITY(klon,zrhol(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:))
	998
	999	DO i = 1, klon
	1000
	1001	IF (rneb(i,k).GT.0.0) THEN
	1002
	1003	! Initialization of zpluie, zice and ztot:
	1004	zpluie=0.
	1005	zice=0.
	1006	ztot=0.
	1007
	1008	IF (zneb(i).EQ.seuil_neb) THEN
	1009	ztot = 0.0
	1010	zice = 0.0
	1011	zpluie= 0.0
	1012	ELSE
	1013	! water quantity to remove: zchau (Sundqvist, 1978)
	1014	! same thing for the ice: zfroi (Zender & Kiehl, 1997)
	1015	IF (ptconv(i,k)) THEN ! if convective point
	1016	zcl=cld_lc_con
	1017	zct=1./cld_tau_con
	1018	zfroi = dtime/REAL(ninter)/zdz(i)zoliq(i)velo(i)*zfice(i)
	1019	ELSE
	1020	zcl=cld_lc_lsc
	1021	zct=1./cld_tau_lsc
	1022	zfroi = dtime/REAL(ninter)/zdz(i)zoliq(i) & ! dqice/dt=1/rhod(rhowiceqice)/dz
	1023	velo(i) zfice(i)
	1024	ENDIF
	1025
	1026	! if vertical heterogeneity is taken into account, we use
	1027	! the "true" volume fraction instead of a modified
	1028	! surface fraction (which is larger and artificially
	1029	! reduces the in-cloud water).
	1030	IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
	1031	zchau = zct dtime/REAL(ninter) zoliq(i) &
	1032	(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl)2)) (1.-zfice(i))
	1033	ELSE
	1034	zchau = zct dtime/REAL(ninter) zoliq(i) &
	1035	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl)2)) (1.-zfice(i)) ! dqliq/dt=-qliq/tau(1-exp(-qcin/clw)*2)
	1036	ENDIF
	1037
	1038	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
	1039	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
	1040	ztot = MAX(zpluie + zice,0.0)
	1041
	1042	ENDIF
	1043
	1044	ztot = MIN(ztot,zoliq(i))
	1045	zice = MIN(zice,ztot)
	1046	zpluie = MIN(zpluie,ztot)
	1047
	1048	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
	1049	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
	1050	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
	1051
	1052	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
	1053	radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
	1054
	1055	ENDIF
	1056
	1057	ENDDO ! i = 1,klon
	1058
	1059	ENDDO ! n = 1,niter
	1060
	1061
	1062	! Caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
	1063	DO i=1,klon
	1064	IF (radliq(i,k) .GT. 0) THEN
	1065	radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
	1066	ENDIF
	1067	ENDDO
	1068
	1069
	1070	! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
	1071	! If T<0C, liquid precip are converted into ice, which leads to an increase in
	1072	! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
	1073	! taken into account through a linearization of the equations and by approximating
	1074	! the liquid precipitation process with a "threshold" process. We assume that
	1075	! condensates are not modified during this operation. Liquid precipitation is
	1076	! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
	1077	! Water vapor increases as well
	1078	! Note that compared to fisrtilp, we always assume iflag_bergeron=2
	1079
	1080
	1081
	1082	DO i = 1, klon
	1083
	1084	IF (rneb(i,k) .GT. 0.0) THEN
	1085
	1086	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
	1087	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
	1088	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
	1089	coef1 = rneb(i,k)RLSTT/zcpzdqsdT_raw(i)
	1090	! Computation of DT if all the liquid precip freezes
	1091	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
	1092	! T should not exceed the freezing point
	1093	! that is Delta > RTT-zt(i)
	1094	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
	1095	zt(i) = zt(i) + DeltaT
	1096	! water vaporization due to temp. increase
	1097	Deltaq = rneb(i,k)zdqsdT_raw(i)DeltaT
	1098	! we add this vaporized water to the total vapor and we remove it from the precip
	1099	zq(i) = zq(i) + Deltaq
	1100	! The three "max" lines herebelow protect from rounding errors
	1101	zcond(i) = max( zcond(i)- Deltaq, 0. )
	1102	! liquid precipitation freezes oe evaporates
	1103	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
	1104	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
	1105	! iced water budget
	1106	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
	1107
	1108	d_ql(i,k) = (1-zfice(i))*zoliq(i)
	1109	d_qi(i,k) = zfice(i)*zoliq(i)
	1110
	1111	IF (iflag_evap_prec.EQ.4) THEN
	1112	zrflcld(i) = zrflcld(i)+zqprecl(i) &
	1113	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1114	ziflcld(i) = ziflcld(i)+ zqpreci(i) &
	1115	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1116	znebprecipcld(i) = rneb(i,k)
	1117	zrfl(i) = zrflcld(i) + zrflclr(i)
	1118	zifl(i) = ziflcld(i) + ziflclr(i)
	1119	ELSE
	1120	zrfl(i) = zrfl(i)+ zqprecl(i) &
	1121	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1122	zifl(i) = zifl(i)+ zqpreci(i) &
	1123	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1124	ENDIF
	1125
	1126	ENDIF ! rneb>0
	1127
	1128	ENDDO
	1129
	1130
	1131
	1132
	1133	! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
	1134	! if iflag_evap_pre=4
	1135	IF (iflag_evap_prec.EQ.4) THEN
	1136
	1137	DO i=1, klon
	1138
	1139
	1140	IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
	1141	znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
	1142	(MAX(znebprecipclr(i),seuil_neb)rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)rain_int_min)))
	1143	ELSE
	1144	znebprecipclr(i)=0.
	1145	ENDIF
	1146
	1147
	1148	IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
	1149	znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
	1150	(MAX(znebprecipcld(i),seuil_neb)rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)rain_int_min)))
	1151	ELSE
	1152	znebprecipcld(i)=0.
	1153	ENDIF
	1154
	1155
	1156	ENDDO
	1157
	1158	ENDIF
	1159
	1160	! --------------------------------
	1161	! End of precipitation formation
	1162	! --------------------------------
	1163
	1164
	1165	! Outputs:
	1166	! Precipitation fluxes at layer interfaces
	1167	! and temperature and vapor tendencies
	1168	DO i = 1, klon
	1169	psfl(i,k)=zifl(i)
	1170	prfl(i,k)=zrfl(i)
	1171	d_q(i,k) = zq(i) - q(i,k)
	1172	d_t(i,k) = zt(i) - t(i,k)
	1173	ENDDO
	1174
	1175
	1176
	1177	! ----------------------------------------------------------------
	1178	! P4> Wet scavenging
	1179	! ----------------------------------------------------------------
	1180
	1181
	1182	!Scavenging through nucleation in the layer
	1183
	1184	DO i = 1,klon
	1185
	1186	IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
	1187	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
	1188	ELSE
	1189	beta(i,k) = 0.
	1190	ENDIF
	1191
	1192	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
	1193	* (paprs(i,k)-paprs(i,k+1))/RG
	1194
	1195
	1196	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
	1197
	1198
	1199	IF (t(i,k) .GE. t_glace_min) THEN
	1200	zalpha_tr = a_tr_sca(3)
	1201	ELSE
	1202	zalpha_tr = a_tr_sca(4)
	1203	ENDIF
	1204
	1205
	1206	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	1207	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
	1208	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
	1209
	1210	! Nucleation with a factor of -1 instead of -0.5
	1211	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
	1212	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
	1213
	1214	ENDIF
	1215
	1216	ENDDO
	1217
	1218
	1219	! Scavenging through impaction in the underlying layer
	1220
	1221	DO kk = k-1, 1, -1
	1222
	1223	DO i = 1, klon
	1224
	1225	IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
	1226
	1227	IF (t(i,kk) .GE. t_glace_min) THEN
	1228	zalpha_tr = a_tr_sca(1)
	1229	ELSE
	1230	zalpha_tr = a_tr_sca(2)
	1231	ENDIF
	1232
	1233	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	1234	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
	1235	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
	1236
	1237	ENDIF
	1238
	1239	ENDDO
	1240
	1241	ENDDO
	1242
	1243
	1244	END DO
	1245
	1246	!======================================================================
	1247	! END OF VERTICAL LOOP
	1248	!======================================================================
	1249
	1250
	1251	! Rain or snow at the surface (depending on the first layer temperature)
	1252	DO i = 1, klon
	1253	snow(i) = zifl(i)
	1254	rain(i) = zrfl(i)
	1255	ENDDO
	1256
	1257
	1258
	1259	IF (ncoreczq>0) THEN
	1260	WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
	1261	ENDIF
	1262
	1263
	1264
	1265
	1266	END SUBROUTINE LSCP
	1267	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1268
	1269	END MODULE LSCP_MOD

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

Download in other formats:

Original Format