Context Navigation

lscp_mod.F90 @ 5453

Last change on this file since 5453 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
File size: 50.0 KB

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

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

Download in other formats:

Original Format