Context Navigation

fisrtilp.F90 @ 3871

Last change on this file since 3871 was 3851, checked in by dcugnet, 4 years ago
Update the branch to the current trunk.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 57.6 KB

Rev	Line
[524]	1	!
[1403]	2	! $Id: fisrtilp.F90 3851 2021-02-22 11:44:07Z acozic $
[524]	3	!
[1472]	4	!
	5	SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
[2086]	6	d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow, &
[1742]	7	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
	8	frac_impa, frac_nucl, beta, &
	9	prfl, psfl, rhcl, zqta, fraca, &
[2236]	10	ztv, zpspsk, ztla, zthl, iflag_cld_th, &
[1849]	11	iflag_ice_thermo)
[524]	12
[1472]	13	!
	14	USE dimphy
[2109]	15	USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
[2311]	16	USE print_control_mod, ONLY: prt_level, lunout
[2686]	17	USE cloudth_mod
[2703]	18	USE ioipsl_getin_p_mod, ONLY : getin_p
[2807]	19	USE phys_local_var_mod, ONLY: ql_seri,qs_seri
[2945]	20	USE phys_local_var_mod, ONLY: rneblsvol
[2807]	21	! flag to include modifications to ensure energy conservation (if flag >0)
	22	USE add_phys_tend_mod, only : fl_cor_ebil
[1472]	23	IMPLICIT none
	24	!======================================================================
	25	! Auteur(s): Z.X. Li (LMD/CNRS)
	26	! Date: le 20 mars 1995
	27	! Objet: condensation et precipitation stratiforme.
	28	! schema de nuage
[2500]	29	! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
	30	! et ilp (il pleut, L. Li)
	31	! Principales parties:
[2807]	32	! P0> Thermalisation des precipitations venant de la couche du dessus
[2500]	33	! P1> Evaporation de la precipitation (qui vient du niveau k+1)
	34	! P2> Formation du nuage (en k)
[2807]	35	! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
	36	! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
	37	! les valeurs de T et Q initiales
	38	! P2.A.2> Prise en compte du couplage entre eau condensee et T.
	39	! P2.A.3> Calcul des valeures finales associees a la formation des nuages
	40	! P2.B> Nuage "tout ou rien"
	41	! P2.C> Prise en compte de la Chaleur latente apres formation nuage
[2500]	42	! P3> Formation de la precipitation (en k)
[1472]	43	!======================================================================
[2807]	44	! JLD:
	45	! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
	46	! * fl_cor_ebil doit etre > 0 ;
	47	! fl_cor_ebil= 0 pour reproduire anciens bugs
[1472]	48	!======================================================================
	49	include "YOMCST.h"
	50	include "fisrtilp.h"
[2006]	51	include "nuage.h" ! JBM (3/14)
[1506]	52
[1472]	53	!
[2500]	54	! Principaux inputs:
[1472]	55	!
[2814]	56	REAL, INTENT(IN) :: dtime ! intervalle du temps (s)
	57	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche
	58	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression au milieu de couche
	59	REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
	60	REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! humidite specifique (kg/kg)
	61	LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! points ou le schema de conv. prof. est actif
	62	INTEGER, INTENT(IN) :: iflag_cld_th
	63	INTEGER, INTENT(IN) :: iflag_ice_thermo
[2500]	64	!
[2814]	65	! Inputs lies aux thermiques
	66	!
	67	REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv
	68	REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta, fraca
	69	REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk, ztla
	70	REAL, DIMENSION(klon,klev), INTENT(IN) :: zthl
	71	!
	72	! Input/output
	73	REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! determine la largeur de distribution de vapeur
	74	!
[2500]	75	! Principaux outputs:
	76	!
[2814]	77	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! incrementation de la temperature (K)
	78	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! incrementation de la vapeur d'eau
	79	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! incrementation de l'eau liquide
	80	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! incrementation de l'eau glace
	81	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! fraction nuageuse
	82	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radliq ! eau liquide utilisee dans rayonnements
	83	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! humidite relative en ciel clair
	84	REAL, DIMENSION(klon), INTENT(OUT) :: rain
	85	REAL, DIMENSION(klon), INTENT(OUT) :: snow
	86	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl
	87	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl
[1403]	88
[1472]	89	!AA
	90	! Coeffients de fraction lessivee : pour OFF-LINE
	91	!
[2814]	92	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_nucl
	93	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_1nucl
	94	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_impa
[1472]	95	!
	96	! Fraction d'aerosols lessivee par impaction et par nucleation
	97	! POur ON-LINE
	98	!
[2814]	99	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa
	100	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl
[1472]	101	!AA
[2814]	102	! --------------------------------------------------------------------------------
[1472]	103	!
	104	! Options du programme:
	105	!
[2923]	106	REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
	107	!$OMP THREADPRIVATE(seuil_neb)
[524]	108
[2923]	109
[1472]	110	INTEGER ninter ! sous-intervals pour la precipitation
	111	PARAMETER (ninter=5)
[2923]	112	INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie
	113	!$OMP THREADPRIVATE(iflag_evap_prec)
[1472]	114	!
	115	LOGICAL cpartiel ! condensation partielle
	116	PARAMETER (cpartiel=.TRUE.)
	117	REAL t_coup
	118	PARAMETER (t_coup=234.0)
[2814]	119	REAL DDT0
	120	PARAMETER (DDT0=.01)
	121	REAL ztfondue
	122	PARAMETER (ztfondue=278.15)
	123	! --------------------------------------------------------------------------------
[1472]	124	!
	125	! Variables locales:
	126	!
	127	INTEGER i, k, n, kk
[2923]	128	INTEGER,save::itap=0
	129	!$OMP THREADPRIVATE(itap)
	130
[2814]	131	REAL qsl, qsi
	132	real zct ,zcl
	133	INTEGER ncoreczq
	134	REAL ctot(klon,klev)
[2945]	135	REAL ctot_vol(klon,klev)
[1901]	136	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
[2807]	137	REAL zdqsdT_raw(klon)
[1901]	138	REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
[2814]	139
	140	logical lognormale(klon)
	141	logical ice_thermo
[1901]	142	LOGICAL convergence(klon)
[2086]	143	INTEGER n_i(klon), iter
	144	REAL cste
[2814]	145
	146	real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
	147	real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
	148	real erf
	149	REAL qcloud(klon)
[1901]	150
[1849]	151	REAL zrfl(klon), zrfln(klon), zqev, zqevt
	152	REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
	153	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
	154	REAL zoliqp(klon), zoliqi(klon)
[2006]	155	REAL zt(klon)
	156	! JBM (3/14) nexpo is replaced by exposant_glace
	157	! REAL nexpo ! exponentiel pour glace/eau
	158	! INTEGER, PARAMETER :: nexpo=6
	159	INTEGER exposant_glace_old
	160	REAL t_glace_min_old
[1472]	161	REAL zdz(klon),zrho(klon),ztot , zrhol(klon)
[2923]	162	REAL zchau ,zfroi ,zfice(klon),zneb(klon),znebprecip(klon)
[2814]	163	REAL zmelt, zpluie, zice
[2086]	164	REAL dzfice(klon)
[2415]	165	REAL zsolid
[2466]	166	!!!!
	167	! Variables pour Bergeron
[2807]	168	REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
[2466]	169	REAL zqpreci(klon), zqprecl(klon)
[2807]	170	! Variable pour conservation enegie des precipitations
	171	REAL zmqc(klon)
[1472]	172	!
	173	LOGICAL appel1er
	174	SAVE appel1er
	175	!$OMP THREADPRIVATE(appel1er)
	176	!
[2703]	177	! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max
	178	! iflag_oldbug_fisrtilp=1 ajoute le BUG
	179	INTEGER,SAVE :: iflag_oldbug_fisrtilp=0 !=0 sans bug
	180	!$OMP THREADPRIVATE(iflag_oldbug_fisrtilp)
[1472]	181	!---------------------------------------------------------------
	182	!
	183	!AA Variables traceurs:
	184	!AA Provisoire !!! Parametres alpha du lessivage
	185	!AA A priori on a 4 scavenging # possibles
	186	!
	187	REAL a_tr_sca(4)
	188	save a_tr_sca
	189	!$OMP THREADPRIVATE(a_tr_sca)
	190	!
	191	! Variables intermediaires
	192	!
	193	REAL zalpha_tr
	194	REAL zfrac_lessi
	195	REAL zprec_cond(klon)
	196	!AA
[1742]	197	! RomP >>> 15 nov 2012
	198	REAL beta(klon,klev) ! taux de conversion de l'eau cond
	199	! RomP <<<
[2807]	200	REAL zmair(klon), zcpair, zcpeau
[1472]	201	! Pour la conversion eau-neige
	202	REAL zlh_solid(klon), zm_solid
	203	!---------------------------------------------------------------
	204	!
	205	! Fonctions en ligne:
	206	!
[2500]	207	REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
	208	! (Heymsfield & Donner, 1990)
[1472]	209	REAL zzz
[2807]	210
[1472]	211	include "YOETHF.h"
	212	include "FCTTRE.h"
	213	fallvc (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_con
	214	fallvs (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_lsc
	215	!
	216	DATA appel1er /.TRUE./
	217	!ym
[2086]	218	!CR: pour iflag_ice_thermo=2, on active que la convection
	219	! ice_thermo = iflag_ice_thermo .GE. 1
[2923]	220
[2945]	221	itap=itap+1
	222	znebprecip(:)=0.
[2923]	223
[2945]	224	ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
[1472]	225	zdelq=0.0
[2945]	226	ctot_vol(1:klon,1:klev)=0.0
	227	rneblsvol(1:klon,1:klev)=0.0
[524]	228
[1506]	229	if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
[1472]	230	IF (appel1er) THEN
[2703]	231	CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
[2923]	232	CALL getin_p('iflag_evap_prec',iflag_evap_prec)
	233	CALL getin_p('seuil_neb',seuil_neb)
[2703]	234	write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp
[1472]	235	!
[1575]	236	WRITE(lunout,*) 'fisrtilp, ninter:', ninter
[2923]	237	WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
[1575]	238	WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
[2923]	239
[1472]	240	IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
[1575]	241	WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
	242	WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
[1472]	243	! CALL abort
	244	ENDIF
	245	appel1er = .FALSE.
	246	!
	247	!AA initialiation provisoire
	248	a_tr_sca(1) = -0.5
	249	a_tr_sca(2) = -0.5
	250	a_tr_sca(3) = -0.5
	251	a_tr_sca(4) = -0.5
	252	!
	253	!AA Initialisation a 1 des coefs des fractions lessivees
	254	!
	255	!cdir collapse
	256	DO k = 1, klev
	257	DO i = 1, klon
	258	pfrac_nucl(i,k)=1.
	259	pfrac_1nucl(i,k)=1.
	260	pfrac_impa(i,k)=1.
[1742]	261	beta(i,k)=0. !RomP initialisation
[1472]	262	ENDDO
	263	ENDDO
[524]	264
[1472]	265	ENDIF ! test sur appel1er
	266	!
	267	!MAf Initialisation a 0 de zoliq
	268	! DO i = 1, klon
	269	! zoliq(i)=0.
	270	! ENDDO
	271	! Determiner les nuages froids par leur temperature
	272	! nexpo regle la raideur de la transition eau liquide / eau glace.
	273	!
[2086]	274	!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
[2006]	275	IF (iflag_t_glace.EQ.0) THEN
	276	! ztglace = RTT - 15.0
	277	t_glace_min_old = RTT - 15.0
	278	!AJ<
	279	IF (ice_thermo) THEN
	280	! nexpo = 2
	281	exposant_glace_old = 2
	282	ELSE
	283	! nexpo = 6
	284	exposant_glace_old = 6
	285	ENDIF
[2086]	286
[1849]	287	ENDIF
[2006]	288
[1849]	289	!! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
	290	!! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
	291	!>AJ
[1472]	292	!cc nexpo = 1
	293	!
	294	! Initialiser les sorties:
	295	!
	296	!cdir collapse
	297	DO k = 1, klev+1
	298	DO i = 1, klon
	299	prfl(i,k) = 0.0
	300	psfl(i,k) = 0.0
	301	ENDDO
	302	ENDDO
[524]	303
[1472]	304	!cdir collapse
	305	DO k = 1, klev
	306	DO i = 1, klon
	307	d_t(i,k) = 0.0
	308	d_q(i,k) = 0.0
	309	d_ql(i,k) = 0.0
[2086]	310	d_qi(i,k) = 0.0
[1472]	311	rneb(i,k) = 0.0
	312	radliq(i,k) = 0.0
	313	frac_nucl(i,k) = 1.
	314	frac_impa(i,k) = 1.
	315	ENDDO
	316	ENDDO
	317	DO i = 1, klon
	318	rain(i) = 0.0
	319	snow(i) = 0.0
	320	zoliq(i)=0.
	321	! ENDDO
	322	!
	323	! Initialiser le flux de precipitation a zero
	324	!
	325	! DO i = 1, klon
	326	zrfl(i) = 0.0
[1849]	327	zifl(i) = 0.0
[1472]	328	zneb(i) = seuil_neb
	329	ENDDO
	330	!
	331	!
	332	!AA Pour plus de securite
[524]	333
[1472]	334	zalpha_tr = 0.
	335	zfrac_lessi = 0.
[524]	336
[2500]	337	!AA==================================================================
[1472]	338	!
	339	ncoreczq=0
[2500]	340	! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
[1472]	341	!
	342	DO k = klev, 1, -1
	343	!
[2500]	344	!AA===============================================================
[1472]	345	!
[2500]	346	! Initialisation temperature et vapeur
[1472]	347	DO i = 1, klon
	348	zt(i)=t(i,k)
	349	zq(i)=q(i,k)
	350	ENDDO
	351	!
[2807]	352	! ----------------------------------------------------------------
	353	! P0> Thermalisation des precipitations venant de la couche du dessus
	354	! ----------------------------------------------------------------
[1472]	355	! Calculer la varition de temp. de l'air du a la chaleur sensible
[2807]	356	! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
	357	! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
	358	! des differentes transformations thermodynamiques. Cette masse d'eau doit
	359	! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
	360	! de l'enthalpie de la couche avec la temperature
	361	! Variables calculees ou modifiees:
	362	! - zt: temperature de la cocuhe
	363	! - zmqc: masse de precip qui doit etre thermalisee
[1472]	364	!
	365	IF(k.LE.klevm1) THEN
	366	DO i = 1, klon
	367	!IM
[2807]	368	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
	369	! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
[1472]	370	zcpair=RCPD(1.0+RVTMP2zq(i))
	371	zcpeau=RCPD*RVTMP2
[2807]	372	if (fl_cor_ebil .GT. 0) then
	373	! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm
	374	! pour s'assurer que la precip arrivant au sol aura bien la temperature de la
	375	! derniere couche
	376	zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
	377	! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
	378	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zmqc(i)zcpeau + zcpair*zt(i) ) &
	379	/ (zcpair + zmqc(i)*zcpeau)
	380	else ! si on maintient les anciennes erreurs
[1472]	381	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zrfl(i)dtime*zcpeau &
[2807]	382	+ zmair(i)zcpairzt(i) ) &
	383	/ (zmair(i)zcpair + zrfl(i)dtime*zcpeau)
	384	end if
[1472]	385	ENDDO
[2885]	386	ELSE ! IF(k.LE.klevm1)
	387	DO i = 1, klon
	388	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
	389	zmqc(i) = 0.
	390	ENDDO
[2807]	391	ENDIF ! end IF(k.LE.klevm1)
	392	!
[2500]	393	! ----------------------------------------------------------------
[2807]	394	! P1> Calcul de l'evaporation de la precipitation
[2500]	395	! ----------------------------------------------------------------
[2807]	396	! On evapore une partie des precipitations venant de la maille du dessus.
	397	! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
	398	! ce que la fraction de cette couche qui est sous le nuage soit saturee.
	399	! Variables calculees ou modifiees:
	400	! - zrfl et zifl: flux de precip liquide et glace
	401	! - zq, zt: humidite et temperature de la cocuhe
	402	! - zmqc: masse de precip qui doit etre thermalisee
	403	!
[2923]	404	IF (iflag_evap_prec>=1) THEN
[1472]	405	DO i = 1, klon
[2807]	406	! S'il y a des precipitations
[1849]	407	IF (zrfl(i)+zifl(i).GT.0.) THEN
[2500]	408	! Calcul du qsat
[1472]	409	IF (thermcep) THEN
	410	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
	411	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
	412	zqs(i)=MIN(0.5,zqs(i))
	413	zcor=1./(1.-RETV*zqs(i))
	414	zqs(i)=zqs(i)*zcor
	415	ELSE
	416	IF (zt(i) .LT. t_coup) THEN
	417	zqs(i) = qsats(zt(i)) / pplay(i,k)
	418	ELSE
	419	zqs(i) = qsatl(zt(i)) / pplay(i,k)
	420	ENDIF
	421	ENDIF
[1849]	422	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	423	ENDDO
	424	!AJ<
[2807]	425
[1849]	426	IF (.NOT. ice_thermo) THEN
	427	DO i = 1, klon
[2807]	428	! S'il y a des precipitations
[1849]	429	IF (zrfl(i)+zifl(i).GT.0.) THEN
[2500]	430	! Evap max pour ne pas saturer la fraction sous le nuage
[3493]	431	! Evap max jusqu'à atteindre la saturation dans la partie
[2807]	432	! de la maille qui est sous le nuage de la couche du dessus
	433	!!! On ne tient compte de cette fraction que sous une seule
	434	!!! couche sous le nuage
[1849]	435	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
[2807]	436	! Ajout de la prise en compte des precip a thermiser
	437	! avec petite reecriture
	438	if (fl_cor_ebil .GT. 0) then ! nouveau
[2500]	439	! Calcul de l'evaporation du flux de precip herite
	440	! d'au-dessus
[1849]	441	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
[2807]	442	* zmair(i)/pplay(i,k)zt(i)RD
	443	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
	444
	445	! Seuil pour ne pas saturer la fraction sous le nuage
	446	zqev = MIN (zqev, zqevt)
	447	! Nouveau flux de precip
	448	zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
	449	! Aucun flux liquide pour T < t_coup, on reevapore tout.
	450	IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
	451	zrfln(i)=0.
	452	zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
	453	END IF
	454	! Nouvelle vapeur
	455	zq(i) = zq(i) + zqev
	456	zmqc(i) = zmqc(i)-zqev
	457	! Nouvelle temperature (chaleur latente)
	458	zt(i) = zt(i) - zqev &
	459	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
[2814]	460	!!JLD debut de partie a supprimer a terme
[2807]	461	else ! if (fl_cor_ebil .GT. 0)
	462	! Calcul de l'evaporation du flux de precip herite
	463	! d'au-dessus
	464	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
[1849]	465	* (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
	466	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
	467	* RG*dtime/(paprs(i,k)-paprs(i,k+1))
[2500]	468	! Seuil pour ne pas saturer la fraction sous le nuage
[1849]	469	zqev = MIN (zqev, zqevt)
[2500]	470	! Nouveau flux de precip
[1849]	471	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
	472	/RG/dtime
[2500]	473	! Aucun flux liquide pour T < t_coup
[1849]	474	IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
[2500]	475	! Nouvelle vapeur
[1849]	476	zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
	477	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[2500]	478	! Nouvelle temperature (chaleur latente)
[1849]	479	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	480	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	481	* RLVTT/RCPD/(1.0+RVTMP2*zq(i))
[2807]	482	end if ! if (fl_cor_ebil .GT. 0)
[2814]	483	!!JLD fin de partie a supprimer a terme
[1849]	484	zrfl(i) = zrfln(i)
	485	zifl(i) = 0.
	486	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	487	ENDDO
	488	!
	489	ELSE ! (.NOT. ice_thermo)
[2807]	490	! ================================
	491	! Avec thermodynamique de la glace
	492	! ================================
[1849]	493	DO i = 1, klon
	494	!AJ<
[2807]	495	! S'il y a des precipitations
	496	IF (zrfl(i)+zifl(i).GT.0.) THEN
[2923]	497
	498	IF (iflag_evap_prec==1) THEN
	499	znebprecip(i)=zneb(i)
	500	ELSE
	501	znebprecip(i)=MAX(zneb(i),znebprecip(i))
	502	ENDIF
[1849]	503
[2807]	504	! Evap max pour ne pas saturer la fraction sous le nuage
[2923]	505	zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
[524]	506
[2807]	507	!JAM
	508	! On differencie qsat pour l'eau et la glace
	509	! Si zdelta=1. --> glace
	510	! Si zdelta=0. --> eau liquide
[2500]	511
	512	! Calcul du qsat par rapport a l'eau liquide
[1849]	513	qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
	514	qsl= MIN(0.5,qsl)
	515	zcor= 1./(1.-RETV*qsl)
	516	qsl= qsl*zcor
	517
[2807]	518	! Calcul de l'evaporation du flux de precip venant du dessus
[2500]	519	! Formulation en racine du flux de precip
	520	! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
[2962]	521	IF (iflag_evap_prec==3) THEN
	522	zqevt = znebprecip(i)coef_eva(1.0-zq(i)/qsl) &
	523	*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
	524	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	525	ELSE
[1849]	526	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
	527	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
[2962]	528	ENDIF
	529
	530
[1849]	531	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
	532	RGdtime/(paprs(i,k)-paprs(i,k+1))
[2500]	533
	534	! Calcul du qsat par rapport a la glace
[1849]	535	qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
	536	qsi= MIN(0.5,qsi)
	537	zcor= 1./(1.-RETV*qsi)
	538	qsi= qsi*zcor
[1472]	539
[2500]	540	! Calcul de la sublimation du flux de precip solide herite
	541	! d'au-dessus
[2962]	542	IF (iflag_evap_prec==3) THEN
	543	zqevti = znebprecip(i)coef_eva(1.0-zq(i)/qsi) &
	544	*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
	545	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	546	ELSE
[1849]	547	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
	548	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
[2962]	549	ENDIF
[1849]	550	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
	551	RGdtime/(paprs(i,k)-paprs(i,k+1))
	552
[2807]	553	!JAM
	554	! Limitation de l'evaporation. On s'assure qu'on ne sature pas
	555	! la fraction de la couche sous le nuage sinon on repartit zqev0
	556	! en conservant la proportion liquide / glace
[1849]	557
	558	IF (zqevt+zqevti.GT.zqev0) THEN
[2807]	559	zqev=zqev0*zqevt/(zqevt+zqevti)
	560	zqevi=zqev0*zqevti/(zqevt+zqevti)
[1849]	561	ELSE
[2807]	562	!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
	563	! liquides et solides meme si on ne sature pas la couche.
	564	! A mon avis, le test est inutile, et il faudrait tout remplacer par:
	565	! zqev=zqevt
	566	! zqevi=zqevti
[1849]	567	IF (zqevt+zqevti.GT.0.) THEN
[2807]	568	zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
	569	zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
[1849]	570	ELSE
	571	zqev=0.
	572	zqevi=0.
	573	ENDIF
	574	ENDIF
[2500]	575	! Nouveaux flux de precip liquide et solide
[1849]	576	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
	577	/RG/dtime)
	578	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
	579	/RG/dtime)
[2500]	580
	581	! Mise a jour de la vapeur, temperature et flux de precip
[1849]	582	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	583	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[2807]	584	if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
	585	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	586	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[1849]	587	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	588	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
[2807]	589	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
	590	+ (zifln(i)-zifl(i)) &
	591	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	592	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
	593	else ! sans correction thermalisation des precips
	594	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	595	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
[1849]	596	* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
	597	+ (zifln(i)-zifl(i)) &
	598	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	599	* RLSTT/RCPD/(1.0+RVTMP2*zq(i))
[2807]	600	end if
	601	! Nouvelles vaeleurs des precips liquides et solides
[1849]	602	zrfl(i) = zrfln(i)
	603	zifl(i) = zifln(i)
[2923]	604	! print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
[1849]	605
[2086]	606	!CR ATTENTION: deplacement de la fonte de la glace
[2466]	607	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
	608	!!! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 !!!!!!!!! jyg
	609	!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
	610	zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! jyg
[2807]	611	! fraction de la precip solide qui est fondue
[2086]	612	zmelt = MIN(MAX(zmelt,0.),1.)
[2500]	613	! Fusion de la glace
[2086]	614	zrfl(i)=zrfl(i)+zmelt*zifl(i)
[2807]	615	if (fl_cor_ebil .LE. 0) then
	616	! the following line should not be here. Indeed, if zifl is modified
	617	! now, zifl(i)*zmelt is no more the amount of ice that has melt
	618	! and therefore the change in temperature computed below is wrong
	619	zifl(i)=zifl(i)*(1.-zmelt)
	620	end if
[2500]	621	! Chaleur latente de fusion
[2807]	622	if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
[2086]	623	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
[2807]	624	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
	625	else ! sans correction thermalisation des precips
	626	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
[2086]	627	RLMLT/RCPD/(1.0+RVTMP2zq(i))
[2807]	628	end if
	629	if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
	630	zifl(i)=zifl(i)*(1.-zmelt)
	631	end if
[2086]	632
[2923]	633	ELSE
	634	! Si on n'a plus de pluies, on reinitialise a 0 la farcion
	635	! sous nuageuse utilisee pour la pluie.
	636	znebprecip(i)=0.
[1849]	637	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
[1472]	638	ENDDO
[1849]	639
	640	ENDIF ! (.NOT. ice_thermo)
	641
[2500]	642	! ----------------------------------------------------------------
	643	! Fin evaporation de la precipitation
	644	! ----------------------------------------------------------------
[2923]	645	ENDIF ! (iflag_evap_prec>=1)
[1472]	646	!
	647	! Calculer Qs et L/Cp*dQs/dT:
	648	!
	649	IF (thermcep) THEN
	650	DO i = 1, klon
[524]	651	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
	652	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
[2807]	653	if (fl_cor_ebil .GT. 0) then ! nouveau
	654	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
	655	else
[524]	656	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
[2807]	657	endif
[524]	658	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
	659	zqs(i) = MIN(0.5,zqs(i))
	660	zcor = 1./(1.-RETV*zqs(i))
	661	zqs(i) = zqs(i)*zcor
	662	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
[2807]	663	zdqsdT_raw(i) = zdqs(i)* &
	664	& RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTTzdelta)
[1472]	665	ENDDO
	666	ELSE
	667	DO i = 1, klon
	668	IF (zt(i).LT.t_coup) THEN
	669	zqs(i) = qsats(zt(i))/pplay(i,k)
	670	zdqs(i) = dqsats(zt(i),zqs(i))
	671	ELSE
	672	zqs(i) = qsatl(zt(i))/pplay(i,k)
	673	zdqs(i) = dqsatl(zt(i),zqs(i))
	674	ENDIF
	675	ENDDO
	676	ENDIF
	677	!
	678	! Determiner la condensation partielle et calculer la quantite
	679	! de l'eau condensee:
	680	!
[1901]	681	!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
[2086]	682	! if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
	683	! write(,) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
	684	! " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
	685	! stop
	686	! endif
[1403]	687
[2500]	688	! ----------------------------------------------------------------
	689	! P2> Formation du nuage
	690	! ----------------------------------------------------------------
[2807]	691	! Variables calculees:
	692	! rneb : fraction nuageuse
	693	! zcond : eau condensee moyenne dans la maille.
	694	! rhcl: humidite relative ciel-clair
	695	! zt : temperature de la maille
	696	! ----------------------------------------------------------------
	697	!
[1472]	698	IF (cpartiel) THEN
[2807]	699	! -------------------------
	700	! P2.A> Nuage fractionnaire
	701	! -------------------------
[1472]	702	!
	703	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
	704	! nuageuse a partir des PDF de Sandrine Bony.
	705	! rneb : fraction nuageuse
	706	! zqn : eau totale dans le nuage
	707	! zcond : eau condensee moyenne dans la maille.
[3493]	708	! on prend en compte le réchauffement qui diminue la partie
[1472]	709	! condensee
	710	!
	711	! Version avec les raqts
[524]	712
[2807]	713	! ----------------------------------------------------------------
	714	! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
	715	! ----------------------------------------------------------------
[1472]	716	if (iflag_pdf.eq.0) then
[524]	717
[2500]	718	! version creneau de (Li, 1998)
[524]	719	do i=1,klon
[1472]	720	zdelq = min(ratqs(i,k),0.99) * zq(i)
	721	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
	722	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
[524]	723	enddo
	724
[2807]	725	else ! if (iflag_pdf.eq.0)
	726	! ----------------------------------------------------------------
	727	! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
	728	! les valeurs de T et Q initiales
	729	! ----------------------------------------------------------------
[524]	730	do i=1,klon
	731	if(zq(i).lt.1.e-15) then
[1472]	732	ncoreczq=ncoreczq+1
	733	zq(i)=1.e-15
[524]	734	endif
[1472]	735	enddo
[1403]	736
[2236]	737	if (iflag_cld_th>=5) then
[1403]	738
[2696]	739	if (iflag_cloudth_vert<=2) then
[2686]	740	call cloudth(klon,klev,k,ztv, &
[1472]	741	zq,zqta,fraca, &
[3493]	742	qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
[1472]	743	ratqs,zqs,t)
[3493]	744	elseif (iflag_cloudth_vert>=3 .and. iflag_cloudth_vert<=5) then
[2686]	745	call cloudth_v3(klon,klev,k,ztv, &
	746	zq,zqta,fraca, &
[3493]	747	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
[2686]	748	ratqs,zqs,t)
[3493]	749	!----------------------------------
	750	!Version these Jean Jouhaud, Decembre 2018
	751	!----------------------------------
	752	elseif (iflag_cloudth_vert==6) then
	753	call cloudth_v6(klon,klev,k,ztv, &
	754	zq,zqta,fraca, &
	755	qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
	756	ratqs,zqs,t)
	757
[2686]	758	endif
[1472]	759	do i=1,klon
[1403]	760	rneb(i,k)=ctot(i,k)
[2945]	761	rneblsvol(i,k)=ctot_vol(i,k)
[1403]	762	zqn(i)=qcloud(i)
[1472]	763	enddo
[1403]	764
[1472]	765	endif
	766
[2236]	767	if (iflag_cld_th <= 4) then
[1472]	768	lognormale = .true.
[2236]	769	elseif (iflag_cld_th >= 6) then
[1472]	770	! lognormale en l'absence des thermiques
	771	lognormale = fraca(:,k) < 1e-10
	772	else
[3493]	773	! Dans le cas iflag_cld_th=5, on prend systématiquement la
[1472]	774	! bi-gaussienne
	775	lognormale = .false.
	776	end if
	777
[2500]	778	!CR: variation de qsat avec T en presence de glace ou non
[2086]	779	!initialisations
[1472]	780	do i=1,klon
[2086]	781	DT(i) = 0.
	782	n_i(i)=0
[1901]	783	Tbef(i)=zt(i)
[2086]	784	qlbef(i)=0.
	785	enddo
	786
[2807]	787	! ----------------------------------------------------------------
	788	! P2.A.2> Prise en compte du couplage entre eau condensee et T.
	789	! Calcul des grandeurs nuageuses en tenant compte de l'effet de
	790	! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
	791	! qui en dependent. Ce changement de temperature est provisoire, et
	792	! la valeur definitive sera calcule plus tard.
	793	! Variables calculees:
	794	! rneb : nebulosite
	795	! zcond: eau condensee en moyenne dans la maille
	796	! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
	797	! pas clair, il n'y a probablement pas de prise en compte de l'effet de
	798	! T sur qsat
	799	! ----------------------------------------------------------------
[2086]	800
	801	!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
	802	if (iflag_fisrtilp_qsat.ge.0) then
[2500]	803	! Iteration pour condensation avec variation de qsat(T)
	804	! -----------------------------------------------------
[2086]	805	do iter=1,iflag_fisrtilp_qsat+1
	806
	807	do i=1,klon
	808	! do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
[2807]	809	! !! convergence = .true. tant que l'on n'a pas converge !!
	810	! ------------------------------
[2086]	811	convergence(i)=abs(DT(i)).gt.DDT0
	812	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
[2807]	813	! si on n'a pas converge
	814	!
	815	! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
	816	! ---------------------------------------------------------------
	817	! Variables calculees:
	818	! rneb : nebulosite
[3851]	819	! zqn : eau totale dans le nuage (in cloud total water content)
[2807]	820	! zcond: eau condensee en moyenne dans la maille
	821	! rhcl: humidite relative ciel-clair
	822	!
	823	Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
[2086]	824	if (.not.ice_thermo) then
	825	zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
	826	else
	827	if (iflag_t_glace.eq.0) then
	828	zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
[2507]	829	else if (iflag_t_glace.ge.1) then
[2703]	830	if (iflag_oldbug_fisrtilp.EQ.0) then
	831	zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
	832	else
[2807]	833	!avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
[2703]	834	zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
	835	endif
[2086]	836	endif
	837	endif
[2807]	838	! Calcul de rneb, qzn et zcond pour les PDF lognormales
[2086]	839	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
[2807]	840	if (fl_cor_ebil .GT. 0) then
	841	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
	842	else
[2086]	843	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
[2807]	844	end if
[2086]	845	zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
	846	zqs(i) = MIN(0.5,zqs(i))
	847	zcor = 1./(1.-RETV*zqs(i))
	848	zqs(i) = zqs(i)*zcor
	849	zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
[1472]	850	zpdf_sig(i)=ratqs(i,k)*zq(i)
	851	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
	852	zpdf_delta(i)=log(zq(i)/zqs(i))
	853	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
	854	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
	855	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
	856	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
	857	zpdf_e1(i)=1.-erf(zpdf_e1(i))
	858	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
	859	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
	860	zpdf_e2(i)=1.-erf(zpdf_e2(i))
[1901]	861
	862	if (zpdf_e1(i).lt.1.e-10) then
	863	rneb(i,k)=0.
	864	zqn(i)=zqs(i)
	865	else
	866	rneb(i,k)=0.5*zpdf_e1(i)
	867	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
	868	endif
	869
[2956]	870	! If vertical heterogeneity, change fraction by volume as well
	871	if (iflag_cloudth_vert>=3) then
	872	ctot_vol(i,k)=rneb(i,k)
	873	rneblsvol(i,k)=ctot_vol(i,k)
	874	endif
	875
[2086]	876	endif !convergence
[2956]	877
[2086]	878	enddo ! boucle en i
	879
[2807]	880	! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
	881	! due a la condensation.
	882	! ---------------------------------------------------------------
	883	! Variables calculees:
	884	! DT : variation de temperature due a la condensation
	885
[2086]	886	if (.not. ice_thermo) then
[2807]	887	! --------------------------
[2086]	888	do i=1,klon
	889	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
	890
[1901]	891	qlbef(i)=max(0.,zqn(i)-zqs(i))
[2807]	892	if (fl_cor_ebil .GT. 0) then
	893	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
	894	else
[1901]	895	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
[2807]	896	end if
[1901]	897	denom=1.+rneb(i,k)*zdqs(i)
	898	DT(i)=num/denom
[2086]	899	n_i(i)=n_i(i)+1
	900	endif
	901	enddo
[1403]	902
[2807]	903	else ! if (.not. ice_thermo)
	904	! --------------------------
[2507]	905	if (iflag_t_glace.ge.1) then
[2109]	906	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
[1472]	907	endif
[1411]	908
[2086]	909	do i=1,klon
	910	if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
	911
	912	if (iflag_t_glace.eq.0) then
	913	zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
	914	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
	915	zfice(i) = zfice(i)**exposant_glace_old
[2807]	916	dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
	917	& / (t_glace_min_old - RTT)
[1901]	918	endif
[2086]	919
[2969]	920	if (iflag_t_glace.ge.1.and.zfice(i)>0.) then
[2807]	921	dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
	922	& / (t_glace_min - t_glace_max)
[2086]	923	endif
	924
	925	if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
	926	dzfice(i)=0.
	927	endif
[1411]	928
[2086]	929	if (zfice(i).lt.1) then
	930	cste=RLVTT
	931	else
	932	cste=RLSTT
	933	endif
	934
	935	qlbef(i)=max(0.,zqn(i)-zqs(i))
[2807]	936	if (fl_cor_ebil .GT. 0) then
	937	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
	938	& +zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
	939	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
	940	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
	941	& qlbef(i)dzfice(i)
	942	else
	943	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
	944	& +zfice(i)RLSTT)/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
	945	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
[2086]	946	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))rneb(i,k)qlbef(i)dzfice(i)
[2807]	947	end if
[2086]	948	DT(i)=num/denom
	949	n_i(i)=n_i(i)+1
	950
	951	endif ! fin convergence
	952	enddo ! fin boucle i
	953
	954	endif !ice_thermo
	955
[2500]	956	enddo ! iter=1,iflag_fisrtilp_qsat+1
	957	! Fin d'iteration pour condensation avec variation de qsat(T)
	958	! -----------------------------------------------------------
[2807]	959	endif ! if (iflag_fisrtilp_qsat.ge.0)
	960	! ----------------------------------------------------------------
	961	! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
	962	! ----------------------------------------------------------------
[524]	963
	964	endif ! iflag_pdf
	965
[2086]	966	! if (iflag_fisrtilp_qsat.eq.-1) then
[2500]	967	!------------------------------------------
	968	!CR: ATTENTION option fausse mais a existe:
	969	! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
	970	! activer les lignes suivantes:
[2086]	971	IF (1.eq.0) THEN
	972	DO i=1,klon
[1146]	973	IF (rneb(i,k) .LE. 0.0) THEN
	974	zqn(i) = 0.0
	975	rneb(i,k) = 0.0
	976	zcond(i) = 0.0
	977	rhcl(i,k)=zq(i)/zqs(i)
	978	ELSE IF (rneb(i,k) .GE. 1.0) THEN
	979	zqn(i) = zq(i)
[1901]	980	rneb(i,k) = 1.0
	981	zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
[1146]	982	rhcl(i,k)=1.0
	983	ELSE
[1901]	984	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
[1146]	985	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
	986	ENDIF
[2500]	987	ENDDO
	988	ENDIF
	989	!------------------------------------------
[1901]	990
[2086]	991	! ELSE
[2807]	992	! ----------------------------------------------------------------
	993	! P2.A.3> Calcul des valeures finales associees a la formation des nuages
	994	! Variables calculees:
	995	! rneb : nebulosite
	996	! zcond: eau condensee en moyenne dans la maille
	997	! zq : eau vapeur dans la maille
	998	! zt : temperature de la maille
	999	! rhcl: humidite relative ciel-clair
	1000	! ----------------------------------------------------------------
	1001	!
	1002	! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
	1003	! Calcule de l'eau condensee moyenne dans la maille (zcond),
	1004	! et de l'humidite relative ciel-clair (rhcl)
[1901]	1005	DO i=1,klon
	1006	IF (rneb(i,k) .LE. 0.0) THEN
	1007	zqn(i) = 0.0
	1008	rneb(i,k) = 0.0
	1009	zcond(i) = 0.0
	1010	rhcl(i,k)=zq(i)/zqs(i)
	1011	ELSE IF (rneb(i,k) .GE. 1.0) THEN
	1012	zqn(i) = zq(i)
	1013	rneb(i,k) = 1.0
	1014	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
	1015	rhcl(i,k)=1.0
	1016	ELSE
	1017	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
	1018	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
	1019	ENDIF
	1020	ENDDO
[2956]	1021	! If vertical heterogeneity, change fraction by volume as well
	1022	if (iflag_cloudth_vert>=3) then
	1023	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
	1024	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
	1025	endif
[1901]	1026
[2086]	1027	! ENDIF
[1901]	1028
[2500]	1029	ELSE ! de IF (cpartiel)
[2807]	1030	! -------------------------
	1031	! P2.B> Nuage "tout ou rien"
	1032	! -------------------------
	1033	! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
[1472]	1034	DO i = 1, klon
	1035	IF (zq(i).GT.zqs(i)) THEN
	1036	rneb(i,k) = 1.0
	1037	ELSE
	1038	rneb(i,k) = 0.0
	1039	ENDIF
	1040	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
	1041	ENDDO
[2807]	1042	ENDIF ! de IF (cpartiel)
[1472]	1043	!
[2500]	1044	! Mise a jour vapeur d'eau
[2807]	1045	! -------------------------
[1472]	1046	DO i = 1, klon
	1047	zq(i) = zq(i) - zcond(i)
	1048	! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
	1049	ENDDO
[1849]	1050	!AJ<
[2807]	1051	! ------------------------------------
	1052	! P2.C> Prise en compte de la Chaleur latente apres formation nuage
[2500]	1053	! -------------------------------------
[2807]	1054	! Variable calcule:
	1055	! zt : temperature de la maille
	1056	!
[1849]	1057	IF (.NOT. ice_thermo) THEN
[1901]	1058	if (iflag_fisrtilp_qsat.lt.1) then
	1059	DO i = 1, klon
	1060	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
	1061	ENDDO
	1062	else if (iflag_fisrtilp_qsat.gt.0) then
	1063	DO i= 1, klon
[2807]	1064	if (fl_cor_ebil .GT. 0) then
	1065	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
	1066	else
[1901]	1067	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
[2807]	1068	end if
[1901]	1069	ENDDO
	1070	endif
[1849]	1071	ELSE
[2507]	1072	if (iflag_t_glace.ge.1) then
[2109]	1073	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
[1901]	1074	endif
[2006]	1075	if (iflag_fisrtilp_qsat.lt.1) then
	1076	DO i = 1, klon
[2109]	1077	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
[2086]	1078	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
	1079	! t_glace_max, exposant_glace)
	1080	if (iflag_t_glace.eq.0) then
[2223]	1081	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
[2086]	1082	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
	1083	zfice(i) = zfice(i)**exposant_glace_old
	1084	endif
[2006]	1085	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
	1086	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
	1087	ENDDO
	1088	else
	1089	DO i=1, klon
[2109]	1090	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
[2086]	1091	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
	1092	! t_glace_max, exposant_glace)
	1093	if (iflag_t_glace.eq.0) then
[2223]	1094	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
[2086]	1095	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
	1096	zfice(i) = zfice(i)**exposant_glace_old
	1097	endif
[2807]	1098	if (fl_cor_ebil .GT. 0) then
	1099	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
	1100	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
	1101	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
	1102	else
[2006]	1103	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
[2807]	1104	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
	1105	end if
[2006]	1106	ENDDO
	1107	endif
	1108	! print*,zt(i),zrfl(i),zifl(i),'temp1'
	1109	ENDIF
[1849]	1110	!>AJ
[2807]	1111
[2500]	1112	! ----------------------------------------------------------------
	1113	! P3> Formation des precipitations
	1114	! ----------------------------------------------------------------
[1472]	1115	!
	1116	! Partager l'eau condensee en precipitation et eau liquide nuageuse
	1117	!
[2500]	1118
	1119	! Initialisation de zoliq (eau condensee moyenne dans la maille)
[1472]	1120	DO i = 1, klon
	1121	IF (rneb(i,k).GT.0.0) THEN
	1122	zoliq(i) = zcond(i)
	1123	zrho(i) = pplay(i,k) / zt(i) / RD
	1124	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
[1849]	1125	ENDIF
	1126	ENDDO
	1127	!AJ<
	1128	IF (.NOT. ice_thermo) THEN
[2006]	1129	IF (iflag_t_glace.EQ.0) THEN
	1130	DO i = 1, klon
	1131	IF (rneb(i,k).GT.0.0) THEN
	1132	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
	1133	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
	1134	zfice(i) = zfice(i)**exposant_glace_old
	1135	! zfice(i) = zfice(i)**nexpo
	1136	!! zfice(i)=0.
	1137	ENDIF
	1138	ENDDO
	1139	ELSE ! of IF (iflag_t_glace.EQ.0)
[2109]	1140	CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
[2086]	1141	! DO i = 1, klon
	1142	! IF (rneb(i,k).GT.0.0) THEN
[2109]	1143	! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
[2086]	1144	! zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
	1145	! t_glace_max, exposant_glace)
	1146	! ENDIF
	1147	! ENDDO
[2006]	1148	ENDIF
[1849]	1149	ENDIF
[2500]	1150
	1151	! Calcul de radliq (eau condensee pour le rayonnement)
	1152	! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
	1153	! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
	1154	! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
	1155	! transmise au rayonnement;
	1156	! ----------------------------------------------------------------
[1849]	1157	DO i = 1, klon
	1158	IF (rneb(i,k).GT.0.0) THEN
[1472]	1159	zneb(i) = MAX(rneb(i,k), seuil_neb)
[1849]	1160	! zt(i) = zt(i)+zcond(i)zfice(i)RLMLT/RCPD/(1.0+RVTMP2*zq(i))
	1161	! print*,zt(i),'fractionglace'
	1162	!>AJ
[1472]	1163	radliq(i,k) = zoliq(i)/REAL(ninter+1)
	1164	ENDIF
	1165	ENDDO
	1166	!
	1167	DO n = 1, ninter
	1168	DO i = 1, klon
	1169	IF (rneb(i,k).GT.0.0) THEN
	1170	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
[1855]	1171	! Initialization of zpluie and zice:
	1172	zpluie=0
	1173	zice=0
[1472]	1174	IF (zneb(i).EQ.seuil_neb) THEN
	1175	ztot = 0.0
	1176	ELSE
[2500]	1177	! quantite d'eau a eliminer: zchau (Sundqvist, 1978)
	1178	! meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
[1472]	1179	if (ptconv(i,k)) then
	1180	zcl =cld_lc_con
	1181	zct =1./cld_tau_con
	1182	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
	1183	fallvc(zrhol(i)) zfice(i)
	1184	else
	1185	zcl =cld_lc_lsc
	1186	zct =1./cld_tau_lsc
	1187	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
	1188	fallvs(zrhol(i)) zfice(i)
	1189	endif
[2945]	1190
	1191	! si l'heterogeneite verticale est active, on utilise
	1192	! la fraction volumique "vraie" plutot que la fraction
	1193	! surfacique modifiee, qui est plus grande et reduit
	1194	! sinon l'eau in-cloud de facon artificielle
	1195	if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
	1196	zchau = zct dtime/REAL(ninter) zoliq(i) &
	1197	(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl )2)) (1.-zfice(i))
	1198	else
	1199	zchau = zct dtime/REAL(ninter) zoliq(i) &
	1200	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )2)) (1.-zfice(i))
	1201	endif
[1849]	1202	!AJ<
	1203	IF (.NOT. ice_thermo) THEN
	1204	ztot = zchau + zfroi
	1205	ELSE
	1206	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
	1207	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
	1208	ztot = zpluie + zice
	1209	ENDIF
	1210	!>AJ
[1472]	1211	ztot = MAX(ztot ,0.0)
	1212	ENDIF
	1213	ztot = MIN(ztot,zoliq(i))
[1849]	1214	!AJ<
	1215	! zoliqp = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
	1216	! zoliqi = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
[2807]	1217	!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
	1218	! temporaires et ne doivent pas etre calcule (alors qu'elles le sont
	1219	! si iflag_bergeron <> 2
	1220	! A SUPPRIMER A TERME
[1849]	1221	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
	1222	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
[1472]	1223	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
[1849]	1224	!>AJ
[1472]	1225	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
	1226	ENDIF
[2466]	1227	ENDDO ! i = 1,klon
	1228	ENDDO ! n = 1,ninter
[2807]	1229
[2500]	1230	! ----------------------------------------------------------------
[1472]	1231	!
[2466]	1232	IF (.NOT. ice_thermo) THEN
[1849]	1233	DO i = 1, klon
	1234	IF (rneb(i,k).GT.0.0) THEN
[1472]	1235	d_ql(i,k) = zoliq(i)
	1236	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
	1237	* (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
[1849]	1238	ENDIF
	1239	ENDDO
	1240	ELSE
[2466]	1241	!
	1242	!CR&JYG<
	1243	! On prend en compte l'effet Bergeron dans les flux de precipitation :
	1244	! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
	1245	! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
	1246	! et les precipitations est grossierement pris en compte en linearisant les equations
	1247	! et en approximant le processus de precipitation liquide par un processus a seuil.
[3493]	1248	! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
[2466]	1249	! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
	1250	! Le condensat precipitant solide est augmente.
	1251	! La vapeur d'eau est augmentee.
	1252	!
	1253	IF ((iflag_bergeron .EQ. 2)) THEN
	1254	DO i = 1, klon
	1255	IF (rneb(i,k) .GT. 0.0) THEN
	1256	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
	1257	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
[2807]	1258	if (fl_cor_ebil .GT. 0) then
	1259	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
	1260	coef1 = rneb(i,k)RLSTT/zcpzdqsdT_raw(i)
	1261	! Calcul de DT si toute les precips liquides congelent
	1262	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
	1263	! On ne veut pas que T devienne superieur a la temp. de congelation.
	1264	! donc que Delta > RTT-zt(i
	1265	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
	1266	zt(i) = zt(i) + DeltaT
	1267	! Eau vaporisee du fait de l'augmentation de T
	1268	Deltaq = rneb(i,k)zdqsdT_raw(i)DeltaT
	1269	! on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
	1270	zq(i) = zq(i) + Deltaq
	1271	! Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
	1272	zcond(i) = max( zcond(i)- Deltaq, 0. )
	1273	! precip liquide qui congele ou qui s'evapore
	1274	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
	1275	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
	1276	! bilan eau glacee
	1277	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
	1278	else ! if (fl_cor_ebil .GT. 0)
	1279	! ancien calcul
[2466]	1280	zcp=RCPD(1.0+RVTMP2(zq(i)+zcond(i)))
	1281	coef1 = RLMLT*zdqs(i)/RLVTT
	1282	DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
	1283	zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
	1284	zqprecl(i) = max( zqprecl(i) - zcp/RLMLT(1.+coef1)DeltaT, 0. )
	1285	zcond(i) = max( zcond(i) - zcp/RLVTTzdqs(i)DeltaT, 0. )
	1286	zq(i) = zq(i) + zcp/RLVTTzdqs(i)DeltaT
	1287	zt(i) = zt(i) + DeltaT
[2807]	1288	end if ! if (fl_cor_ebil .GT. 0)
[2466]	1289	ENDIF ! rneb(i,k) .GT. 0.0
	1290	ENDDO
	1291	DO i = 1, klon
	1292	IF (rneb(i,k).GT.0.0) THEN
	1293	d_ql(i,k) = (1-zfice(i))*zoliq(i)
	1294	d_qi(i,k) = zfice(i)*zoliq(i)
	1295	zrfl(i) = zrfl(i)+ zqprecl(i) &
	1296	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1297	zifl(i) = zifl(i)+ zqpreci(i) &
	1298	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1299	ENDIF
	1300	ENDDO
	1301	!!
	1302	ELSE ! iflag_bergeron
	1303	!>CR&JYG
	1304	!!
[1849]	1305	DO i = 1, klon
	1306	IF (rneb(i,k).GT.0.0) THEN
[2086]	1307	!CR on prend en compte la phase glace
[2807]	1308	!JLD inutile car on ne passe jamais ici si .not.ice_thermo
	1309	! if (.not.ice_thermo) then
	1310	! d_ql(i,k) = zoliq(i)
	1311	! d_qi(i,k) = 0.
	1312	! else
[2086]	1313	d_ql(i,k) = (1-zfice(i))*zoliq(i)
	1314	d_qi(i,k) = zfice(i)*zoliq(i)
[2807]	1315	! endif
[1849]	1316	!AJ<
	1317	zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
	1318	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1319	zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
	1320	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1321	! zrfl(i) = zrfl(i)+ zpluie &
	1322	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1323	! zifl(i) = zifl(i)+ zice &
	1324	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1325
[2415]	1326	!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
[2466]	1327	IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
[2415]	1328	zsolid = zrfl(i)
	1329	zifl(i) = zifl(i)+zrfl(i)
	1330	zrfl(i) = 0.
[2807]	1331	if (fl_cor_ebil .GT. 0) then
[2415]	1332	zt(i)=zt(i)+zsolid(RGdtime)/(paprs(i,k)-paprs(i,k+1)) &
[2807]	1333	(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
	1334	else
	1335	zt(i)=zt(i)+zsolid(RGdtime)/(paprs(i,k)-paprs(i,k+1)) &
[2415]	1336	(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2zq(i))
[2807]	1337	end if
[2466]	1338	ENDIF ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
[2415]	1339	!RC
	1340
[2466]	1341	ENDIF ! rneb(i,k).GT.0.0
[1849]	1342	ENDDO
	1343
[2466]	1344	ENDIF ! iflag_bergeron .EQ. 2
	1345	ENDIF ! .NOT. ice_thermo
	1346
[2086]	1347	!CR: la fonte est faite au debut
	1348	! IF (ice_thermo) THEN
	1349	! DO i = 1, klon
	1350	! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
	1351	! zmelt = MIN(MAX(zmelt,0.),1.)
	1352	! zrfl(i)=zrfl(i)+zmelt*zifl(i)
	1353	! zifl(i)=zifl(i)*(1.-zmelt)
[1849]	1354	! print*,zt(i),'octavio1'
[2086]	1355	! zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
	1356	! RLMLT/RCPD/(1.0+RVTMP2zq(i))
[1849]	1357	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
[2086]	1358	! ENDDO
	1359	! ENDIF
[1849]	1360
	1361
	1362	IF (.NOT. ice_thermo) THEN
	1363	DO i = 1, klon
	1364	IF (zt(i).LT.RTT) THEN
[1472]	1365	psfl(i,k)=zrfl(i)
[1849]	1366	ELSE
[1472]	1367	prfl(i,k)=zrfl(i)
[1849]	1368	ENDIF
	1369	ENDDO
	1370	ELSE
	1371	! JAM*************************************************
[2500]	1372	! Revoir partie ci-dessous: a quoi servent psfl et prfl?
[1849]	1373	! *****************************************************
	1374
	1375	DO i = 1, klon
	1376	! IF (zt(i).LT.RTT) THEN
	1377	psfl(i,k)=zifl(i)
	1378	! ELSE
	1379	prfl(i,k)=zrfl(i)
	1380	! ENDIF
	1381	!>AJ
	1382	ENDDO
	1383	ENDIF
[2500]	1384	! ----------------------------------------------------------------
	1385	! Fin de formation des precipitations
	1386	! ----------------------------------------------------------------
[1472]	1387	!
	1388	! Calculer les tendances de q et de t:
	1389	!
	1390	DO i = 1, klon
	1391	d_q(i,k) = zq(i) - q(i,k)
	1392	d_t(i,k) = zt(i) - t(i,k)
	1393	ENDDO
	1394	!
	1395	!AA--------------- Calcul du lessivage stratiforme -------------
[524]	1396
[1472]	1397	DO i = 1,klon
	1398	!
[1742]	1399	if(zcond(i).gt.zoliq(i)+1.e-10) then
	1400	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
	1401	else
	1402	beta(i,k) = 0.
	1403	endif
[1472]	1404	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
	1405	* (paprs(i,k)-paprs(i,k+1))/RG
	1406	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
	1407	!AA lessivage nucleation LMD5 dans la couche elle-meme
[2006]	1408	IF (iflag_t_glace.EQ.0) THEN
	1409	if (t(i,k) .GE. t_glace_min_old) THEN
[1472]	1410	zalpha_tr = a_tr_sca(3)
	1411	else
	1412	zalpha_tr = a_tr_sca(4)
	1413	endif
[2006]	1414	ELSE ! of IF (iflag_t_glace.EQ.0)
	1415	if (t(i,k) .GE. t_glace_min) THEN
	1416	zalpha_tr = a_tr_sca(3)
	1417	else
	1418	zalpha_tr = a_tr_sca(4)
	1419	endif
	1420	ENDIF
[1472]	1421	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	1422	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
	1423	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
	1424	!
	1425	! nucleation avec un facteur -1 au lieu de -0.5
	1426	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
	1427	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
	1428	ENDIF
	1429	!
	1430	ENDDO ! boucle sur i
	1431	!
	1432	!AA Lessivage par impaction dans les couches en-dessous
	1433	DO kk = k-1, 1, -1
[524]	1434	DO i = 1, klon
[1472]	1435	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
[2006]	1436	IF (iflag_t_glace.EQ.0) THEN
	1437	if (t(i,kk) .GE. t_glace_min_old) THEN
[1472]	1438	zalpha_tr = a_tr_sca(1)
	1439	else
	1440	zalpha_tr = a_tr_sca(2)
	1441	endif
[2006]	1442	ELSE ! of IF (iflag_t_glace.EQ.0)
	1443	if (t(i,kk) .GE. t_glace_min) THEN
	1444	zalpha_tr = a_tr_sca(1)
	1445	else
	1446	zalpha_tr = a_tr_sca(2)
	1447	endif
	1448	ENDIF
[1472]	1449	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	1450	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
	1451	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
	1452	ENDIF
[524]	1453	ENDDO
[1472]	1454	ENDDO
	1455	!
[2500]	1456	!AA===============================================================
	1457	! FIN DE LA BOUCLE VERTICALE
[1472]	1458	end DO
	1459	!
[2500]	1460	!AA==================================================================
[1472]	1461	!
	1462	! Pluie ou neige au sol selon la temperature de la 1ere couche
	1463	!
[2086]	1464	!CR: si la thermo de la glace est active, on calcule zifl directement
	1465	IF (.NOT.ice_thermo) THEN
[1472]	1466	DO i = 1, klon
	1467	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
[1849]	1468	!AJ<
[2086]	1469	! snow(i) = zrfl(i)
[1849]	1470	snow(i) = zrfl(i)+zifl(i)
	1471	!>AJ
[1472]	1472	zlh_solid(i) = RLSTT-RLVTT
	1473	ELSE
	1474	rain(i) = zrfl(i)
	1475	zlh_solid(i) = 0.
	1476	ENDIF
	1477	ENDDO
[2086]	1478
	1479	ELSE
	1480	DO i = 1, klon
	1481	snow(i) = zifl(i)
	1482	rain(i) = zrfl(i)
	1483	ENDDO
	1484
	1485	ENDIF
[1472]	1486	!
	1487	! For energy conservation : when snow is present, the solification
	1488	! latent heat is considered.
[2086]	1489	!CR: si thermo de la glace, neige deja prise en compte
	1490	IF (.not.ice_thermo) THEN
[1472]	1491	DO k = 1, klev
	1492	DO i = 1, klon
	1493	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
[2807]	1494	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
[1472]	1495	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
[2807]	1496	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair(i))
[1472]	1497	END DO
	1498	END DO
[2086]	1499	ENDIF
[1472]	1500	!
[883]	1501
[1472]	1502	if (ncoreczq>0) then
[1575]	1503	WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
[1472]	1504	endif
	1505
	1506	END SUBROUTINE fisrtilp

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

Download in other formats:

Original Format