Context Navigation

conlmd.F90 @ 2284

Last change on this file since 2284 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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: 64.1 KB

Rev	Line
[1992]	1
[524]	2	! $Header$
[1992]	3
	4	SUBROUTINE conlmd(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, rain, snow, &
	5	ibas, itop)
	6	USE dimphy
	7	IMPLICIT NONE
	8	! ======================================================================
	9	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	10	! Objet: Schema de convection utilis'e dans le modele du LMD
	11	! Ajustement humide (Manabe) + Ajustement convectif (Kuo)
	12	! ======================================================================
	13	! ym#include "dimensions.h"
	14	! ym#include "dimphy.h"
	15	include "YOMCST.h"
	16	include "YOETHF.h"
	17
	18	! Arguments:
	19
	20	REAL dtime ! pas d'integration (s)
	21	REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
	22	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	23	REAL t(klon, klev) ! temperature (K)
	24	REAL q(klon, klev) ! humidite specifique (kg/kg)
	25	REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)
	26
	27	REAL d_t(klon, klev) ! incrementation temperature
	28	REAL d_q(klon, klev) ! incrementation humidite
	29	REAL rain(klon) ! pluies (mm/s)
	30	REAL snow(klon) ! neige (mm/s)
	31	INTEGER ibas(klon) ! niveau du bas
	32	INTEGER itop(klon) ! niveau du haut
	33
	34	LOGICAL usekuo ! utiliser convection profonde (schema Kuo)
	35	PARAMETER (usekuo=.TRUE.)
	36
	37	REAL d_t_bis(klon, klev)
	38	REAL d_q_bis(klon, klev)
	39	REAL rain_bis(klon)
	40	REAL snow_bis(klon)
	41	INTEGER ibas_bis(klon)
	42	INTEGER itop_bis(klon)
	43	REAL d_ql(klon, klev), d_ql_bis(klon, klev)
	44	REAL rneb(klon, klev), rneb_bis(klon, klev)
	45
	46	INTEGER i, k
	47	REAL zlvdcp, zlsdcp, zdelta, zz, za, zb
	48
	49	! cc CALL fiajh ! ancienne version de Convection Manabe
	50	CALL conman & ! nouvelle version de Convection
	51	! Manabe
	52	(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, ibas, itop)
	53
	54	IF (usekuo) THEN
	55	! cc CALL fiajc ! ancienne version de Convection Kuo
	56	CALL conkuo & ! nouvelle version de Convection
	57	! Kuo
	58	(dtime, paprs, pplay, t, q, conv_q, d_t_bis, d_q_bis, d_ql_bis, &
	59	rneb_bis, rain_bis, snow_bis, ibas_bis, itop_bis)
	60	DO k = 1, klev
[524]	61	DO i = 1, klon
[1992]	62	d_t(i, k) = d_t(i, k) + d_t_bis(i, k)
	63	d_q(i, k) = d_q(i, k) + d_q_bis(i, k)
	64	d_ql(i, k) = d_ql(i, k) + d_ql_bis(i, k)
	65	END DO
	66	END DO
	67	DO i = 1, klon
	68	rain(i) = rain(i) + rain_bis(i)
	69	snow(i) = snow(i) + snow_bis(i)
	70	ibas(i) = min(ibas(i), ibas_bis(i))
	71	itop(i) = max(itop(i), itop_bis(i))
	72	END DO
	73	END IF
	74
	75	! L'eau liquide convective est dispersee dans l'air:
	76
	77	DO k = 1, klev
	78	DO i = 1, klon
	79	zlvdcp = rlvtt/rcpd/(1.0+rvtmp2*q(i,k))
	80	zlsdcp = rlstt/rcpd/(1.0+rvtmp2*q(i,k))
	81	zdelta = max(0., sign(1.,rtt-t(i,k)))
	82	zz = d_ql(i, k) ! re-evap. de l'eau liquide
	83	zb = max(0.0, zz)
	84	za = -max(0.0, zz)(zlvdcp(1.-zdelta)+zlsdcp*zdelta)
	85	d_t(i, k) = d_t(i, k) + za
	86	d_q(i, k) = d_q(i, k) + zb
	87	END DO
	88	END DO
	89
	90	RETURN
	91	END SUBROUTINE conlmd
	92	SUBROUTINE conman(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, &
	93	snow, ibas, itop)
	94	USE dimphy
	95	IMPLICIT NONE
	96	! ======================================================================
	97	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19970324
	98	! Objet: ajustement humide convectif avec la possibilite de faire
	99	! l'ajustement sur une fraction de la maille.
	100	! Methode: On impose une distribution uniforme pour la vapeur d'eau
	101	! au sein d'une maille. On applique la procedure d'ajustement
	102	! successivement a la totalite, 75%, 50%, 25% et 5% de la maille
	103	! jusqu'a ce que l'ajustement a lieu. J'espere que ceci augmente
	104	! les activites convectives et corrige le biais "trop froid et sec"
	105	! du modele.
	106	! ======================================================================
	107	! ym#include "dimensions.h"
	108	! ym#include "dimphy.h"
	109	include "YOMCST.h"
	110
	111	REAL dtime ! pas d'integration (s)
	112	REAL t(klon, klev) ! temperature (K)
	113	REAL q(klon, klev) ! humidite specifique (kg/kg)
	114	REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
	115	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	116
	117	REAL d_t(klon, klev) ! incrementation temperature
	118	REAL d_q(klon, klev) ! incrementation humidite
	119	REAL d_ql(klon, klev) ! incrementation eau liquide
	120	REAL rneb(klon, klev) ! nebulosite
	121	REAL rain(klon) ! pluies (mm/s)
	122	REAL snow(klon) ! neige (mm/s)
	123	INTEGER ibas(klon) ! niveau du bas
	124	INTEGER itop(klon) ! niveau du haut
	125
	126	LOGICAL afaire(klon) ! .TRUE. implique l'ajustement
	127	LOGICAL accompli(klon) ! .TRUE. si l'ajustement est effectif
	128
	129	INTEGER nb ! nombre de sous-fractions a considere
	130	PARAMETER (nb=1)
	131	! cc PARAMETER (nb=3)
	132
	133	REAL ratqs ! largeur de la distribution pour vapeur d'eau
	134	PARAMETER (ratqs=0.05)
	135
	136	REAL w_q(klon, klev)
	137	REAL w_d_t(klon, klev), w_d_q(klon, klev), w_d_ql(klon, klev)
	138	REAL w_rneb(klon, klev)
	139	REAL w_rain(klon), w_snow(klon)
	140	INTEGER w_ibas(klon), w_itop(klon)
	141	REAL zq1, zq2
	142	INTEGER i, k, n
	143
	144	REAL t_coup
	145	PARAMETER (t_coup=234.0)
	146	REAL zdp1, zdp2
	147	REAL zqs1, zqs2, zdqs1, zdqs2
	148	REAL zgamdz
	149	REAL zflo ! flotabilite
	150	REAL zsat ! sur-saturation
	151	REAL zdelta, zcor, zcvm5
	152	LOGICAL imprim
	153
	154	INTEGER ncpt
	155	SAVE ncpt
	156	!$OMP THREADPRIVATE(ncpt)
	157	REAL frac(nb) ! valeur de la maille fractionnelle
	158	SAVE frac
	159	!$OMP THREADPRIVATE(frac)
	160	INTEGER opt_cld(nb) ! option pour le modele nuageux
	161	SAVE opt_cld
	162	!$OMP THREADPRIVATE(opt_cld)
	163	LOGICAL appel1er
	164	SAVE appel1er
	165	!$OMP THREADPRIVATE(appel1er)
	166
	167	! Fonctions thermodynamiques:
	168
	169	include "YOETHF.h"
	170	include "FCTTRE.h"
	171
	172	DATA frac/1.0/
	173	DATA opt_cld/4/
	174	! cc DATA frac / 1.0, 0.50, 0.25/
	175	! cc DATA opt_cld / 4, 4, 4/
	176
	177	DATA appel1er/.TRUE./
	178	DATA ncpt/0/
	179
	180	IF (appel1er) THEN
	181	PRINT *, 'conman, nb:', nb
	182	PRINT *, 'conman, frac:', frac
	183	PRINT *, 'conman, opt_cld:', opt_cld
	184	appel1er = .FALSE.
	185	END IF
	186
	187	! Initialiser les sorties a zero:
	188
	189	DO k = 1, klev
	190	DO i = 1, klon
	191	d_t(i, k) = 0.0
	192	d_q(i, k) = 0.0
	193	d_ql(i, k) = 0.0
	194	rneb(i, k) = 0.0
	195	END DO
	196	END DO
	197	DO i = 1, klon
	198	ibas(i) = klev
	199	itop(i) = 1
	200	rain(i) = 0.0
	201	snow(i) = 0.0
	202	END DO
	203
	204	! S'il n'y a pas d'instabilite conditionnelle,
	205	! pas la penne de se fatiguer:
	206
	207	DO i = 1, klon
	208	afaire(i) = .FALSE.
	209	END DO
	210	DO k = 1, klev - 1
	211	DO i = 1, klon
	212	IF (thermcep) THEN
	213	zdelta = max(0., sign(1.,rtt-t(i,k)))
	214	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	215	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
	216	zqs1 = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
	217	zqs1 = min(0.5, zqs1)
	218	zcor = 1./(1.-retv*zqs1)
	219	zqs1 = zqs1*zcor
	220	zdqs1 = foede(t(i,k), zdelta, zcvm5, zqs1, zcor)
	221
	222	zdelta = max(0., sign(1.,rtt-t(i,k+1)))
	223	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	224	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k+1))
	225	zqs2 = r2es*foeew(t(i,k+1), zdelta)/pplay(i, k+1)
	226	zqs2 = min(0.5, zqs2)
	227	zcor = 1./(1.-retv*zqs2)
	228	zqs2 = zqs2*zcor
	229	zdqs2 = foede(t(i,k+1), zdelta, zcvm5, zqs2, zcor)
	230	ELSE
	231	IF (t(i,k)<t_coup) THEN
	232	zqs1 = qsats(t(i,k))/pplay(i, k)
	233	zdqs1 = dqsats(t(i,k), zqs1)
	234
	235	zqs2 = qsats(t(i,k+1))/pplay(i, k+1)
	236	zdqs2 = dqsats(t(i,k+1), zqs2)
	237	ELSE
	238	zqs1 = qsatl(t(i,k))/pplay(i, k)
	239	zdqs1 = dqsatl(t(i,k), zqs1)
	240
	241	zqs2 = qsatl(t(i,k+1))/pplay(i, k+1)
	242	zdqs2 = dqsatl(t(i,k+1), zqs2)
	243	END IF
	244	END IF
	245	zdp1 = paprs(i, k) - paprs(i, k+1)
	246	zdp2 = paprs(i, k+1) - paprs(i, k+2)
	247	zgamdz = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd(rd(t(i, &
	248	k)zdp1+t(i,k+1)zdp2)/(zdp1+zdp2)+rlvtt(zqs1zdp1+zqs2*zdp2)/(zdp1+ &
	249	zdp2))/(1.0+(zdqs1zdp1+zdqs2zdp2)/(zdp1+zdp2))
	250	zflo = t(i, k) + zgamdz - t(i, k+1)
	251	zsat = (q(i,k)-zqs1)zdp1 + (q(i,k+1)-zqs2)zdp2
	252	IF (zflo>0.0) afaire(i) = .TRUE.
	253	! erreur IF (zflo.GT.0.0 .AND. zsat.GT.0.0) afaire(i) = .TRUE.
	254	END DO
	255	END DO
	256
	257	imprim = mod(ncpt, 48) == 0
	258	DO n = 1, nb
	259
	260	DO k = 1, klev
[524]	261	DO i = 1, klon
[1992]	262	IF (afaire(i)) THEN
	263	zq1 = q(i, k)*(1.0-ratqs)
	264	zq2 = q(i, k)*(1.0+ratqs)
	265	w_q(i, k) = zq2 - frac(n)/2.0*(zq2-zq1)
	266	END IF
	267	END DO
	268	END DO
	269
	270	CALL conmanv(dtime, paprs, pplay, t, w_q, afaire, opt_cld(n), w_d_t, &
	271	w_d_q, w_d_ql, w_rneb, w_rain, w_snow, w_ibas, w_itop, accompli, &
	272	imprim)
	273	DO k = 1, klev
[524]	274	DO i = 1, klon
[1992]	275	IF (afaire(i) .AND. accompli(i)) THEN
	276	d_t(i, k) = w_d_t(i, k)*frac(n)
	277	d_q(i, k) = w_d_q(i, k)*frac(n)
	278	d_ql(i, k) = w_d_ql(i, k)*frac(n)
	279	IF (nint(w_rneb(i,k))==1) rneb(i, k) = frac(n)
	280	END IF
	281	END DO
	282	END DO
	283	DO i = 1, klon
[524]	284	IF (afaire(i) .AND. accompli(i)) THEN
[1992]	285	rain(i) = w_rain(i)*frac(n)
	286	snow(i) = w_snow(i)*frac(n)
	287	ibas(i) = min(ibas(i), w_ibas(i))
	288	itop(i) = max(itop(i), w_itop(i))
	289	END IF
	290	END DO
	291	DO i = 1, klon
	292	IF (afaire(i) .AND. accompli(i)) afaire(i) = .FALSE.
	293	END DO
	294
	295	END DO
	296
	297	ncpt = ncpt + 1
	298
	299	RETURN
	300	END SUBROUTINE conman
	301	SUBROUTINE conmanv(dtime, paprs, pplay, t, q, afaire, opt_cld, d_t, d_q, &
	302	d_ql, rneb, rain, snow, ibas, itop, accompli, imprim)
	303	USE dimphy
	304	IMPLICIT NONE
	305	! ======================================================================
	306	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	307	! Objet: ajustement humide (convection proposee par Manabe).
	308	! Pour une colonne verticale, il peut avoir plusieurs blocs
	309	! necessitant l'ajustement. ibas est le bas du plus bas bloc
	310	! et itop est le haut du plus haut bloc
	311	! ======================================================================
	312	! ym#include "dimensions.h"
	313	! ym#include "dimphy.h"
	314	include "YOMCST.h"
	315
	316	! Arguments:
	317
	318	REAL dtime ! pas d'integration (s)
	319	REAL t(klon, klev) ! temperature (K)
	320	REAL q(klon, klev) ! humidite specifique (kg/kg)
	321	REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
	322	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	323	INTEGER opt_cld ! comment traiter l'eau liquide
	324	LOGICAL afaire(klon) ! .TRUE. si le point est a faire (Input)
	325	LOGICAL imprim ! .T. pour imprimer quelques diagnostiques
	326
	327	REAL d_t(klon, klev) ! incrementation temperature
	328	REAL d_q(klon, klev) ! incrementation humidite
	329	REAL d_ql(klon, klev) ! incrementation eau liquide
	330	REAL rneb(klon, klev) ! nebulosite
	331	REAL rain(klon) ! pluies (mm/s)
	332	REAL snow(klon) ! neige (mm/s)
	333	INTEGER ibas(klon) ! niveau du bas
	334	INTEGER itop(klon) ! niveau du haut
	335	LOGICAL accompli(klon) ! .TRUE. si l'ajustement a eu lieu (Output)
	336
	337	! Quelques options:
	338
	339	LOGICAL new_top ! re-calculer sommet quand re-ajustement est fait
	340	PARAMETER (new_top=.FALSE.)
	341	LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
	342	PARAMETER (evap_prec=.TRUE.)
	343	REAL coef_eva
	344	PARAMETER (coef_eva=1.0E-05)
	345	REAL t_coup
	346	PARAMETER (t_coup=234.0)
	347	REAL seuil_vap
	348	PARAMETER (seuil_vap=1.0E-10)
	349	LOGICAL old_tau ! implique precip nulle, si vrai.
	350	PARAMETER (old_tau=.FALSE.)
	351	REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
	352	REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
	353	PARAMETER (dpmin=0.15, tomax=0.97)
	354	REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
	355	PARAMETER (dpmax=0.30, tomin=0.05)
	356	REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
	357	PARAMETER (deep_sig=0.50, deep_to=0.05)
	358	LOGICAL exigent ! implique un calcul supplementaire pour Qs
	359	PARAMETER (exigent=.FALSE.)
	360
	361	INTEGER kbase
	362	PARAMETER (kbase=0)
	363
	364	! Variables locales:
	365
	366	INTEGER nexpo
	367	INTEGER i, k, k1min, k1max, k2min, k2max, is
	368	REAL zgamdz(klon, klev-1)
	369	REAL zt(klon, klev), zq(klon, klev)
	370	REAL zqs(klon, klev), zdqs(klon, klev)
	371	REAL zqmqsdp(klon, klev)
	372	REAL ztnew(klon, klev), zqnew(klon, klev)
	373	REAL zcond(klon), zvapo(klon), zrapp(klon)
	374	REAL zrfl(klon), zrfln, zqev, zqevt
	375	REAL zsat(klon) ! sur-saturation
	376	REAL zflo(klon) ! flotabilite
	377	REAL za(klon), zb(klon), zc(klon)
	378	INTEGER k1(klon), k2(klon)
	379	REAL zdelta, zcor, zcvm5
	380	REAL delp(klon, klev)
	381	LOGICAL possible(klon), todo(klon), etendre(klon)
	382	LOGICAL aller(klon), todobis(klon)
	383	REAL zalfa
	384	INTEGER nbtodo, nbdone
	385
	386	! Fonctions thermodynamiques:
	387
	388	include "YOETHF.h"
	389	include "FCTTRE.h"
	390
	391	DO k = 1, klev
	392	DO i = 1, klon
	393	delp(i, k) = paprs(i, k) - paprs(i, k+1)
	394	END DO
	395	END DO
	396
	397	! Initialiser les sorties a zero
	398
	399	DO k = 1, klev
	400	DO i = 1, klon
	401	d_t(i, k) = 0.0
	402	d_q(i, k) = 0.0
	403	d_ql(i, k) = 0.0
	404	rneb(i, k) = 0.0
	405	END DO
	406	END DO
	407	DO i = 1, klon
	408	ibas(i) = klev
	409	itop(i) = 1
	410	rain(i) = 0.0
	411	snow(i) = 0.0
	412	accompli(i) = .FALSE.
	413	END DO
	414
	415	! Preparations
	416
	417	DO k = 1, klev
	418	DO i = 1, klon
[524]	419	IF (afaire(i)) THEN
[1992]	420	zt(i, k) = t(i, k)
	421	zq(i, k) = q(i, k)
	422
	423	! Calculer Qs et L/Cp*dQs/dT
	424
	425	IF (thermcep) THEN
	426	zdelta = max(0., sign(1.,rtt-zt(i,k)))
	427	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	428	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i,k))
	429	zqs(i, k) = r2es*foeew(zt(i,k), zdelta)/pplay(i, k)
	430	zqs(i, k) = min(0.5, zqs(i,k))
	431	zcor = 1./(1.-retv*zqs(i,k))
	432	zqs(i, k) = zqs(i, k)*zcor
	433	zdqs(i, k) = foede(zt(i,k), zdelta, zcvm5, zqs(i,k), zcor)
	434	ELSE
	435	IF (zt(i,k)<t_coup) THEN
	436	zqs(i, k) = qsats(zt(i,k))/pplay(i, k)
	437	zdqs(i, k) = dqsats(zt(i,k), zqs(i,k))
	438	ELSE
	439	zqs(i, k) = qsatl(zt(i,k))/pplay(i, k)
	440	zdqs(i, k) = dqsatl(zt(i,k), zqs(i,k))
	441	END IF
	442	END IF
	443
	444	! Calculer (q-qs)*dp
	445	zqmqsdp(i, k) = (zq(i,k)-zqs(i,k))*delp(i, k)
	446	END IF
	447	END DO
	448	END DO
	449
	450	! -----zgama is the moist convective lapse rate (-dT/dz).
	451	! -----zgamdz(*,k) est la difference minimale autorisee des temperatures
	452	! -----entre deux couches (k et k+1), c.a.d. si T(k+1)-T(k) est inferieur
	453	! -----a zgamdz(*,k), alors ces 2 couches sont instables conditionnellement
	454
	455	DO k = 1, klev - 1
	456	DO i = 1, klon
[524]	457	IF (afaire(i)) THEN
[1992]	458	zgamdz(i, k) = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd(rd(zt( &
	459	i,k)delp(i,k)+zt(i,k+1)delp(i,k+1))/(delp(i,k)+delp(i, &
	460	k+1))+rlvtt(zqs(i,k)delp(i,k)+zqs(i,k+1)*delp(i,k+1))/(delp(i, &
	461	k)+delp(i,k+1)))/(1.0+(zdqs(i,k)delp(i,k)+zdqs(i,k+1)delp(i, &
	462	k+1))/(delp(i,k)+delp(i,k+1)))
	463	END IF
	464	END DO
	465	END DO
	466
	467	! On cherche la presence simultanee d'instabilite conditionnelle
	468	! et de sur-saturation. Sinon, pas la penne de se fatiguer:
	469
	470	DO i = 1, klon
	471	possible(i) = .FALSE.
	472	END DO
	473	DO k = 2, klev
	474	DO i = 1, klon
[524]	475	IF (afaire(i)) THEN
[1992]	476	zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
	477	zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k-1)
	478	IF (zflo(i)>0.0 .AND. zsat(i)>0.0) possible(i) = .TRUE.
	479	END IF
	480	END DO
	481	END DO
	482
	483	DO i = 1, klon
	484	IF (possible(i)) THEN
	485	k1(i) = kbase
	486	k2(i) = k1(i) + 1
	487	END IF
	488	END DO
	489
	490	810 CONTINUE ! chercher le bas de la colonne a ajuster
	491
	492	k2min = klev
	493	DO i = 1, klon
	494	todo(i) = .FALSE.
	495	aller(i) = .TRUE.
	496	IF (possible(i)) k2min = min(k2min, k2(i))
	497	END DO
	498	IF (k2min==klev) GO TO 860
	499	DO k = k2min, klev - 1
	500	DO i = 1, klon
	501	IF (possible(i) .AND. k>=k2(i) .AND. aller(i)) THEN
	502	zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k+1)
	503	zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k+1)
	504	IF (zflo(i)>0.0 .AND. zsat(i)>0.0) THEN
	505	k1(i) = k
	506	k2(i) = k + 1
	507	todo(i) = .TRUE.
	508	aller(i) = .FALSE.
	509	END IF
	510	END IF
	511	END DO
	512	END DO
	513	DO i = 1, klon
	514	IF (possible(i) .AND. aller(i)) THEN
	515	todo(i) = .FALSE.
	516	k1(i) = klev
	517	k2(i) = klev
	518	END IF
	519	END DO
	520
	521	! CC DO i = 1, klon
	522	! CC IF (possible(i)) THEN
	523	! CC 811 k2(i) = k2(i) + 1
	524	! CC IF (k2(i) .GT. klev) THEN
	525	! CC todo(i) = .FALSE.
	526	! CC GOTO 812
	527	! CC ENDIF
	528	! CC k = k2(i)
	529	! CC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
	530	! CC zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1)
	531	! CC IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 811
	532	! CC k1(i) = k2(i) - 1
	533	! CC todo(i) = .TRUE.
	534	! CC ENDIF
	535	! CC 812 CONTINUE
	536	! CC ENDDO
	537
	538	820 CONTINUE ! chercher le haut de la colonne
	539
	540	k2min = klev
	541	DO i = 1, klon
	542	aller(i) = .TRUE.
	543	IF (todo(i)) k2min = min(k2min, k2(i))
	544	END DO
	545	IF (k2min<klev) THEN
	546	DO k = k2min, klev
[524]	547	DO i = 1, klon
[1992]	548	IF (todo(i) .AND. k>k2(i) .AND. aller(i)) THEN
	549	zsat(i) = zsat(i) + zqmqsdp(i, k)
	550	zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
	551	IF (zflo(i)<=0.0 .OR. zsat(i)<=0.0) THEN
[524]	552	aller(i) = .FALSE.
[1992]	553	ELSE
	554	k2(i) = k
	555	END IF
	556	END IF
	557	END DO
	558	END DO
	559	! error is = 0
	560	! error DO i = 1, klon
	561	! error IF(todo(i).AND.aller(i)) THEN
	562	! error is = is + 1
	563	! error todo(i) = .FALSE.
	564	! error k2(i) = klev
	565	! error ENDIF
	566	! error ENDDO
	567	! error IF (is.GT.0) THEN
	568	! error PRINT*, "Bizard. je pourrais continuer mais j arrete"
	569	! error CALL abort
	570	! error ENDIF
	571	END IF
	572
	573	! CC DO i = 1, klon
	574	! CC IF (todo(i)) THEN
	575	! CC 821 CONTINUE
	576	! CC IF (k2(i) .EQ. klev) GOTO 822
	577	! CC k = k2(i) + 1
	578	! CC zsat(i) = zsat(i) + zqmqsdp(i,k)
	579	! CC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
	580	! CC IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 822
	581	! CC k2(i) = k
	582	! CC GOTO 821
	583	! CC ENDIF
	584	! CC 822 CONTINUE
	585	! CC ENDDO
	586
	587	830 CONTINUE ! faire l'ajustement en sachant k1 et k2
	588
	589	is = 0
	590	DO i = 1, klon
	591	IF (todo(i)) THEN
	592	IF (k2(i)<=k1(i)) is = is + 1
	593	END IF
	594	END DO
	595	IF (is>0) THEN
	596	PRINT *, 'Impossible: k1 trop grand ou k2 trop petit'
	597	PRINT *, 'is=', is
	598	CALL abort
	599	END IF
	600
	601	k1min = klev
	602	k1max = 1
	603	k2max = 1
	604	DO i = 1, klon
	605	IF (todo(i)) THEN
	606	k1min = min(k1min, k1(i))
	607	k1max = max(k1max, k1(i))
	608	k2max = max(k2max, k2(i))
	609	END IF
	610	END DO
	611
	612	DO i = 1, klon
	613	IF (todo(i)) THEN
[524]	614	k = k1(i)
	615	za(i) = 0.
[1992]	616	zb(i) = (rcpd(1.+zdqs(i,k))(zt(i,k)-za(i))-rlvtt*(zqs(i,k)-zq(i, &
	617	k)))*delp(i, k)
	618	zc(i) = delp(i, k)rcpd(1.+zdqs(i,k))
	619	END IF
	620	END DO
	621
	622	DO k = k1min, k2max
	623	DO i = 1, klon
	624	IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) THEN
	625	za(i) = za(i) + zgamdz(i, k-1)
	626	zb(i) = zb(i) + (rcpd(1.+zdqs(i,k))(zt(i,k)-za(i))-rlvtt*(zqs(i, &
	627	k)-zq(i,k)))*delp(i, k)
	628	zc(i) = zc(i) + delp(i, k)rcpd(1.+zdqs(i,k))
	629	END IF
	630	END DO
	631	END DO
	632
	633	DO i = 1, klon
	634	IF (todo(i)) THEN
	635	k = k1(i)
	636	ztnew(i, k) = zb(i)/zc(i)
	637	zqnew(i, k) = zqs(i, k) + (ztnew(i,k)-zt(i,k))rcpd/rlvttzdqs(i, k)
	638	END IF
	639	END DO
	640
	641	DO k = k1min, k2max
	642	DO i = 1, klon
	643	IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) THEN
	644	ztnew(i, k) = ztnew(i, k-1) + zgamdz(i, k-1)
	645	zqnew(i, k) = zqs(i, k) + (ztnew(i,k)-zt(i,k))rcpd/rlvttzdqs(i, k)
	646	END IF
	647	END DO
	648	END DO
	649
	650	! Quantite de condensation produite pendant l'ajustement:
	651
	652	DO i = 1, klon
	653	zcond(i) = 0.0
	654	END DO
	655	DO k = k1min, k2max
	656	DO i = 1, klon
	657	IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
	658	rneb(i, k) = 1.0
	659	zcond(i) = zcond(i) + (zq(i,k)-zqnew(i,k))*delp(i, k)/rg
	660	END IF
	661	END DO
	662	END DO
	663
	664	! Si condensation negative, effort completement perdu:
	665
	666	DO i = 1, klon
	667	IF (todo(i) .AND. zcond(i)<=0.) todo(i) = .FALSE.
	668	END DO
	669
	670	! L'ajustement a ete accompli, meme les calculs accessoires
	671	! ne sont pas encore faits:
	672
	673	DO i = 1, klon
	674	IF (todo(i)) accompli(i) = .TRUE.
	675	END DO
	676
	677	! =====
	678	! Une fois que la condensation a lieu, on doit construire un
	679	! "modele nuageux" pour partager la condensation entre l'eau
	680	! liquide nuageuse et la precipitation (leur rapport toliq
	681	! est calcule selon l'epaisseur nuageuse). Je suppose que
	682	! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
	683	! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
	684	! lineaire entre dpmin et dpmax).
	685	! =====
	686	DO i = 1, klon
	687	IF (todo(i)) THEN
	688	toliq(i) = tomax - ((paprs(i,k1(i))-paprs(i,k2(i)+1))/paprs(i,1)-dpmin) &
	689	*(tomax-tomin)/(dpmax-dpmin)
	690	toliq(i) = max(tomin, min(tomax,toliq(i)))
	691	IF (pplay(i,k2(i))/paprs(i,1)<=deep_sig) toliq(i) = deep_to
	692	IF (old_tau) toliq(i) = 1.0
	693	END IF
	694	END DO
	695	! =====
	696	! On doit aussi determiner la distribution verticale de
	697	! l'eau nuageuse. Plusieurs options sont proposees:
	698
	699	! (0) La condensation precipite integralement (toliq ne sera
	700	! pas utilise).
	701	! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
	702	! a la vapeur d'eau locale.
	703	! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
	704	! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
	705	! est effectivement diminuee pendant le processus d'ajustement.
	706	! (4) Elle est en fonction (lineaire ou exponentielle) de la
	707	! distance (epaisseur en pression) avec le niveau k1 (la couche
	708	! k1 n'aura donc pas d'eau liquide).
	709	! =====
	710
	711	IF (opt_cld==0) THEN
	712
	713	DO i = 1, klon
	714	IF (todo(i)) zrfl(i) = zcond(i)/dtime
	715	END DO
	716
	717	ELSE IF (opt_cld==1) THEN
	718
	719	DO i = 1, klon
	720	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
	721	END DO
	722	DO k = k1min, k2max
[524]	723	DO i = 1, klon
[1992]	724	IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
	725	zqnew(i, k)*delp(i, k)/rg
	726	END DO
	727	END DO
	728	DO i = 1, klon
[524]	729	IF (todo(i)) THEN
[1992]	730	zrapp(i) = toliq(i)*zcond(i)/zvapo(i)
	731	zrapp(i) = max(0., min(1.,zrapp(i)))
	732	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	733	END IF
	734	END DO
	735	DO k = k1min, k2max
[524]	736	DO i = 1, klon
[1992]	737	IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
	738	d_ql(i, k) = d_ql(i, k) + zrapp(i)*zqnew(i, k)
	739	END IF
	740	END DO
	741	END DO
	742
	743	ELSE IF (opt_cld==2) THEN
	744
	745	DO i = 1, klon
	746	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
	747	END DO
	748	DO k = k1min, k2max
[524]	749	DO i = 1, klon
[1992]	750	IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
	751	delp(i, k)/rg
	752	END DO
	753	END DO
	754	DO k = k1min, k2max
[524]	755	DO i = 1, klon
[1992]	756	IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
	757	d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)
	758	END IF
	759	END DO
	760	END DO
	761	DO i = 1, klon
	762	IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	763	END DO
	764
	765	ELSE IF (opt_cld==3) THEN
	766
	767	DO i = 1, klon
	768	IF (todo(i)) zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
	769	END DO
	770	DO k = k1min, k2max
[524]	771	DO i = 1, klon
[1992]	772	IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
	773	max(0.0, zq(i,k)-zqnew(i,k))*delp(i, k)/rg
	774	END DO
	775	END DO
	776	DO k = k1min, k2max
[524]	777	DO i = 1, klon
[1992]	778	IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i) .AND. zvapo(i)>0.0) d_ql(i, &
	779	k) = d_ql(i, k) + toliq(i)zcond(i)/zvapo(i)max(0.0, zq(i,k)-zqnew &
	780	(i,k))
	781	END DO
	782	END DO
	783	DO i = 1, klon
	784	IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	785	END DO
	786
	787	ELSE IF (opt_cld==4) THEN
	788
	789	nexpo = 3
	790	! cc nexpo = 1 ! distribution lineaire
	791
	792	DO i = 1, klon
	793	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
	794	END DO ! (avec ponderation)
	795	DO k = k1min, k2max
[524]	796	DO i = 1, klon
[1992]	797	IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
	798	delp(i, k)/rg(pplay(i,k1(i))-pplay(i,k))*nexpo
	799	END DO
	800	END DO
	801	DO k = k1min, k2max
[524]	802	DO i = 1, klon
[1992]	803	IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) d_ql(i, k) = d_ql(i, &
	804	k) + toliq(i)zcond(i)/zvapo(i)(pplay(i,k1(i))-pplay(i,k))**nexpo
	805	END DO
	806	END DO
	807	DO i = 1, klon
	808	IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	809	END DO
	810
	811	ELSE ! valeur non-prevue pour opt_cld
	812
	813	PRINT *, 'opt_cld est faux:', opt_cld
	814	CALL abort
	815
	816	END IF ! fin de opt_cld
	817
	818	! L'eau precipitante peut etre evaporee:
	819
	820	zalfa = 0.05
	821	IF (evap_prec .AND. (k1max>=2)) THEN
	822	DO k = k1max - 1, 1, -1
[524]	823	DO i = 1, klon
[1992]	824	IF (todo(i) .AND. k<k1(i) .AND. zrfl(i)>0.0) THEN
	825	zqev = max(0.0, (zqs(i,k)-zq(i,k))*zalfa)
	826	zqevt = coef_eva(1.0-zq(i,k)/zqs(i,k))sqrt(zrfl(i))*delp(i, k)/ &
	827	pplay(i, k)zt(i, k)rd/rg
	828	zqevt = max(0.0, min(zqevt,zrfl(i)))rgdtime/delp(i, k)
	829	zqev = min(zqev, zqevt)
	830	zrfln = zrfl(i) - zqev*(delp(i,k))/rg/dtime
	831	zq(i, k) = zq(i, k) - (zrfln-zrfl(i))(rg/(delp(i,k)))dtime
	832	zt(i, k) = zt(i, k) + (zrfln-zrfl(i))*(rg/(delp(i, &
	833	k)))dtimerlvtt/rcpd/(1.0+rvtmp2*zq(i,k))
	834	zrfl(i) = zrfln
	835	END IF
	836	END DO
	837	END DO
	838	END IF
	839
	840	! La temperature de la premiere couche determine la pluie ou la neige:
	841
	842	DO i = 1, klon
	843	IF (todo(i)) THEN
	844	IF (zt(i,1)>rtt) THEN
	845	rain(i) = rain(i) + zrfl(i)
[524]	846	ELSE
[1992]	847	snow(i) = snow(i) + zrfl(i)
	848	END IF
	849	END IF
	850	END DO
	851
	852	! Mise a jour de la temperature et de l'humidite
	853
	854	DO k = k1min, k2max
	855	DO i = 1, klon
	856	IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
	857	zt(i, k) = ztnew(i, k)
	858	zq(i, k) = zqnew(i, k)
	859	END IF
	860	END DO
	861	END DO
	862
	863	! Re-calculer certaines variables pour etendre et re-ajuster la colonne
	864
	865	IF (exigent) THEN
	866	DO k = 1, klev
[524]	867	DO i = 1, klon
[1992]	868	IF (todo(i)) THEN
	869	IF (thermcep) THEN
	870	zdelta = max(0., sign(1.,rtt-zt(i,k)))
	871	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	872	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i,k))
	873	zqs(i, k) = r2es*foeew(zt(i,k), zdelta)/pplay(i, k)
	874	zqs(i, k) = min(0.5, zqs(i,k))
	875	zcor = 1./(1.-retv*zqs(i,k))
	876	zqs(i, k) = zqs(i, k)*zcor
	877	zdqs(i, k) = foede(zt(i,k), zdelta, zcvm5, zqs(i,k), zcor)
	878	ELSE
	879	IF (zt(i,k)<t_coup) THEN
	880	zqs(i, k) = qsats(zt(i,k))/pplay(i, k)
	881	zdqs(i, k) = dqsats(zt(i,k), zqs(i,k))
	882	ELSE
	883	zqs(i, k) = qsatl(zt(i,k))/pplay(i, k)
	884	zdqs(i, k) = dqsatl(zt(i,k), zqs(i,k))
	885	END IF
	886	END IF
	887	END IF
	888	END DO
	889	END DO
	890	END IF
	891
	892	IF (exigent) THEN
	893	DO k = 1, klev - 1
[524]	894	DO i = 1, klon
[1992]	895	IF (todo(i)) THEN
	896	zgamdz(i, k) = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd(rd( &
	897	zt(i,k)delp(i,k)+zt(i,k+1)delp(i,k+1))/(delp(i,k)+delp(i, &
	898	k+1))+rlvtt(zqs(i,k)delp(i,k)+zqs(i,k+1)*delp(i,k+1))/(delp(i, &
	899	k)+delp(i,k+1)))/(1.0+(zdqs(i,k)delp(i,k)+zdqs(i,k+1)delp(i, &
	900	k+1))/(delp(i,k)+delp(i,k+1)))
	901	END IF
	902	END DO
	903	END DO
	904	END IF
	905
	906	! Puisque l'humidite a ete modifiee, on re-fait (q-qs)*dp
	907
	908	DO k = 1, klev
	909	DO i = 1, klon
[524]	910	IF (todo(i)) THEN
[1992]	911	zqmqsdp(i, k) = (zq(i,k)-zqs(i,k))*delp(i, k)
	912	END IF
	913	END DO
	914	END DO
	915
	916	! Verifier si l'on peut etendre le bas de la colonne
	917
	918	DO i = 1, klon
	919	etendre(i) = .FALSE.
	920	END DO
	921
	922	k1max = 1
	923	DO i = 1, klon
	924	IF (todo(i) .AND. k1(i)>(kbase+1)) THEN
	925	k = k1(i)
	926	zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
	927	zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k-1)
	928	! sc voici l'ancienne ligne:
	929	! sc IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) THEN
	930	! sc sylvain: il faut RESPECTER les 2 criteres:
	931	IF (zflo(i)>0.0 .AND. zsat(i)>0.0) THEN
	932	etendre(i) = .TRUE.
	933	k1(i) = k1(i) - 1
	934	k1max = max(k1max, k1(i))
	935	aller(i) = .TRUE.
	936	END IF
	937	END IF
	938	END DO
	939
	940	IF (k1max>(kbase+1)) THEN
	941	DO k = k1max, kbase + 1, -1
[524]	942	DO i = 1, klon
[1992]	943	IF (etendre(i) .AND. k<k1(i) .AND. aller(i)) THEN
	944	zsat(i) = zsat(i) + zqmqsdp(i, k)
	945	zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k+1)
	946	IF (zsat(i)<=0.0 .OR. zflo(i)<=0.0) THEN
[524]	947	aller(i) = .FALSE.
[1992]	948	ELSE
[524]	949	k1(i) = k
[1992]	950	END IF
	951	END IF
	952	END DO
	953	END DO
	954	DO i = 1, klon
	955	IF (etendre(i) .AND. aller(i)) THEN
	956	k1(i) = 1
	957	END IF
	958	END DO
	959	END IF
	960
	961	! CC DO i = 1, klon
	962	! CC IF (etendre(i)) THEN
	963	! CC 840 k = k1(i)
	964	! CC IF (k.GT.1) THEN
	965	! CC zsat(i) = zsat(i) + zqmqsdp(i,k-1)
	966	! CC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
	967	! CC IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN
	968	! CC k1(i) = k - 1
	969	! CC GOTO 840
	970	! CC ENDIF
	971	! CC ENDIF
	972	! CC ENDIF
	973	! CC ENDDO
	974
	975	DO i = 1, klon
	976	todobis(i) = todo(i)
	977	todo(i) = .FALSE.
	978	END DO
	979	is = 0
	980	DO i = 1, klon
	981	IF (etendre(i)) THEN
	982	todo(i) = .TRUE.
	983	is = is + 1
	984	END IF
	985	END DO
	986	IF (is>0) THEN
	987	IF (new_top) THEN
	988	GO TO 820 ! chercher de nouveau le sommet k2
	989	ELSE
	990	GO TO 830 ! supposer que le sommet est celui deja trouve
	991	END IF
	992	END IF
	993
	994	DO i = 1, klon
	995	possible(i) = .FALSE.
	996	END DO
	997	is = 0
	998	DO i = 1, klon
	999	IF (todobis(i) .AND. k2(i)<klev) THEN
	1000	is = is + 1
	1001	possible(i) = .TRUE.
	1002	END IF
	1003	END DO
	1004	IF (is>0) GO TO 810 !on cherche en haut d'autres blocks
	1005	! a ajuster a partir du sommet de la colonne precedente
	1006
	1007	860 CONTINUE ! Calculer les tendances et diagnostiques
	1008	! cc print*, "Apres 860"
	1009
	1010	DO k = 1, klev
	1011	DO i = 1, klon
[524]	1012	IF (accompli(i)) THEN
[1992]	1013	d_t(i, k) = zt(i, k) - t(i, k)
	1014	zq(i, k) = max(zq(i,k), seuil_vap)
	1015	d_q(i, k) = zq(i, k) - q(i, k)
	1016	END IF
	1017	END DO
	1018	END DO
	1019
	1020	DO i = 1, klon
	1021	IF (accompli(i)) THEN
[524]	1022	DO k = 1, klev
[1992]	1023	IF (rneb(i,k)>0.0) THEN
	1024	ibas(i) = k
	1025	GO TO 807
	1026	END IF
	1027	END DO
	1028	807 CONTINUE
	1029	DO k = klev, 1, -1
	1030	IF (rneb(i,k)>0.0) THEN
	1031	itop(i) = k
	1032	GO TO 808
	1033	END IF
	1034	END DO
	1035	808 CONTINUE
	1036	END IF
	1037	END DO
	1038
	1039	IF (imprim) THEN
	1040	nbtodo = 0
	1041	nbdone = 0
	1042	DO i = 1, klon
	1043	IF (afaire(i)) nbtodo = nbtodo + 1
	1044	IF (accompli(i)) nbdone = nbdone + 1
	1045	END DO
	1046	PRINT *, 'nbTodo, nbDone=', nbtodo, nbdone
	1047	END IF
	1048
	1049	RETURN
	1050	END SUBROUTINE conmanv
	1051	SUBROUTINE conkuo(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
	1052	rain, snow, ibas, itop)
	1053	USE dimphy
	1054	IMPLICIT NONE
	1055	! ======================================================================
	1056	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	1057	! Objet: Schema de convection de type Kuo (1965).
	1058	! Cette version du code peut calculer le niveau de depart
	1059	! N.B. version vectorielle (le 6 oct. 1997)
	1060	! ======================================================================
	1061	! ym#include "dimensions.h"
	1062	! ym#include "dimphy.h"
	1063	include "YOMCST.h"
	1064
	1065	! Arguments:
	1066
	1067	REAL dtime ! intervalle du temps (s)
	1068	REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
	1069	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	1070	REAL t(klon, klev) ! temperature (K)
	1071	REAL q(klon, klev) ! humidite specifique
	1072	REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)
	1073
	1074	REAL d_t(klon, klev) ! incrementation temperature
	1075	REAL d_q(klon, klev) ! incrementation humidite
	1076	REAL d_ql(klon, klev) ! incrementation eau liquide
	1077	REAL rneb(klon, klev) ! nebulosite
	1078	REAL rain(klon) ! pluies (mm/s)
	1079	REAL snow(klon) ! neige (mm/s)
	1080	INTEGER itop(klon) ! niveau du sommet
	1081	INTEGER ibas(klon) ! niveau du bas
	1082
	1083	LOGICAL ldcum(klon) ! convection existe
	1084	LOGICAL todo(klon)
	1085
	1086	! Quelsques options:
	1087
	1088	LOGICAL calcfcl ! calculer le niveau de convection libre
	1089	PARAMETER (calcfcl=.TRUE.)
	1090	INTEGER ldepar ! niveau fixe de convection libre
	1091	PARAMETER (ldepar=4)
	1092	INTEGER opt_cld ! comment traiter l'eau liquide
	1093	PARAMETER (opt_cld=4) ! valeur possible: 0, 1, 2, 3 ou 4
	1094	LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
	1095	PARAMETER (evap_prec=.TRUE.)
	1096	REAL coef_eva
	1097	PARAMETER (coef_eva=1.0E-05)
	1098	LOGICAL new_deh ! nouvelle facon de calculer dH
	1099	PARAMETER (new_deh=.FALSE.)
	1100	REAL t_coup
	1101	PARAMETER (t_coup=234.0)
	1102	LOGICAL old_tau ! implique precipitation nulle
	1103	PARAMETER (old_tau=.FALSE.)
	1104	REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
	1105	REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
	1106	PARAMETER (dpmin=0.15, tomax=0.97)
	1107	REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
	1108	PARAMETER (dpmax=0.30, tomin=0.05)
	1109	REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
	1110	PARAMETER (deep_sig=0.50, deep_to=0.05)
	1111
	1112	! Variables locales:
	1113
	1114	INTEGER nexpo
	1115	LOGICAL nuage(klon)
	1116	INTEGER i, k, kbmin, kbmax, khmax
	1117	REAL ztotal(klon, klev), zdeh(klon, klev)
	1118	REAL zgz(klon, klev)
	1119	REAL zqs(klon, klev)
	1120	REAL zdqs(klon, klev)
	1121	REAL ztemp(klon, klev)
	1122	REAL zpres(klon, klev)
	1123	REAL zconv(klon) ! convergence d'humidite
	1124	REAL zvirt(klon) ! convergence virtuelle d'humidite
	1125	REAL zfrac(klon) ! fraction convective
	1126	INTEGER kb(klon), kh(klon)
	1127
	1128	REAL zcond(klon), zvapo(klon), zrapp(klon)
	1129	REAL zrfl(klon), zrfln, zqev, zqevt
	1130	REAL zdelta, zcvm5, zcor
	1131	REAL zvar
	1132
	1133	LOGICAL appel1er
	1134	SAVE appel1er
	1135	!$OMP THREADPRIVATE(appel1er)
	1136
	1137	! Fonctions thermodynamiques
	1138
	1139	include "YOETHF.h"
	1140	include "FCTTRE.h"
	1141
	1142	DATA appel1er/.TRUE./
	1143
	1144	IF (appel1er) THEN
	1145	PRINT *, 'conkuo, calcfcl:', calcfcl
	1146	IF (.NOT. calcfcl) PRINT *, 'conkuo, ldepar:', ldepar
	1147	PRINT *, 'conkuo, opt_cld:', opt_cld
	1148	PRINT *, 'conkuo, evap_prec:', evap_prec
	1149	PRINT *, 'conkuo, new_deh:', new_deh
	1150	appel1er = .FALSE.
	1151	END IF
	1152
	1153	! Initialiser les sorties a zero
	1154
	1155	DO k = 1, klev
	1156	DO i = 1, klon
	1157	d_q(i, k) = 0.0
	1158	d_t(i, k) = 0.0
	1159	d_ql(i, k) = 0.0
	1160	rneb(i, k) = 0.0
	1161	END DO
	1162	END DO
	1163	DO i = 1, klon
	1164	rain(i) = 0.0
	1165	snow(i) = 0.0
	1166	ibas(i) = 0
	1167	itop(i) = 0
	1168	END DO
	1169
	1170	! Calculer la vapeur d'eau saturante Qs et sa derive L/Cp * dQs/dT
	1171
	1172	DO k = 1, klev
	1173	DO i = 1, klon
	1174	IF (thermcep) THEN
	1175	zdelta = max(0., sign(1.,rtt-t(i,k)))
	1176	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1177	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
	1178	zqs(i, k) = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
	1179	zqs(i, k) = min(0.5, zqs(i,k))
	1180	zcor = 1./(1.-retv*zqs(i,k))
	1181	zqs(i, k) = zqs(i, k)*zcor
	1182	zdqs(i, k) = foede(t(i,k), zdelta, zcvm5, zqs(i,k), zcor)
	1183	ELSE
	1184	IF (t(i,k)<t_coup) THEN
	1185	zqs(i, k) = qsats(t(i,k))/pplay(i, k)
	1186	zdqs(i, k) = dqsats(t(i,k), zqs(i,k))
	1187	ELSE
	1188	zqs(i, k) = qsatl(t(i,k))/pplay(i, k)
	1189	zdqs(i, k) = dqsatl(t(i,k), zqs(i,k))
	1190	END IF
	1191	END IF
	1192	END DO
	1193	END DO
	1194
	1195	! Calculer gz (energie potentielle)
	1196
	1197	DO i = 1, klon
	1198	zgz(i, 1) = rdt(i, 1)/(0.5(paprs(i,1)+pplay(i, &
	1199	1)))*(paprs(i,1)-pplay(i,1))
	1200	END DO
	1201	DO k = 2, klev
	1202	DO i = 1, klon
	1203	zgz(i, k) = zgz(i, k-1) + rd0.5(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i &
	1204	,k-1)-pplay(i,k))
	1205	END DO
	1206	END DO
	1207
	1208	! Calculer l'energie statique humide saturee (CpT + gz + LQs)
	1209
	1210	DO k = 1, klev
	1211	DO i = 1, klon
	1212	ztotal(i, k) = rcpdt(i, k) + rlvttzqs(i, k) + zgz(i, k)
	1213	END DO
	1214	END DO
	1215
	1216	! Determiner le niveau de depart et calculer la difference de
	1217	! l'energie statique humide saturee (ztotal) entre la couche
	1218	! de depart et chaque couche au-dessus.
	1219
	1220	IF (calcfcl) THEN
	1221	DO k = 1, klev
[524]	1222	DO i = 1, klon
[1992]	1223	zpres(i, k) = pplay(i, k)
	1224	ztemp(i, k) = t(i, k)
	1225	END DO
	1226	END DO
	1227	CALL kuofcl(ztemp, q, zgz, zpres, ldcum, kb)
	1228	DO i = 1, klon
	1229	IF (ldcum(i)) THEN
	1230	k = kb(i)
	1231	IF (new_deh) THEN
	1232	zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
	1233	ELSE
	1234	zdeh(i, k) = rcpd(t(i,k-1)-t(i,k)) - rd0.5*(t(i,k-1)+t(i,k))/ &
	1235	paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
	1236	rlvtt*(zqs(i,k-1)-zqs(i,k))
	1237	END IF
	1238	zdeh(i, k) = zdeh(i, k)*0.5
	1239	END IF
	1240	END DO
	1241	DO k = 1, klev
[524]	1242	DO i = 1, klon
[1992]	1243	IF (ldcum(i) .AND. k>=(kb(i)+1)) THEN
	1244	IF (new_deh) THEN
	1245	zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
	1246	ELSE
	1247	zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
	1248	rd0.5(t(i,k-1)+t(i,k))/paprs(i, k)* &
	1249	(pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
	1250	END IF
	1251	END IF
	1252	END DO
	1253	END DO
	1254	ELSE
	1255	DO i = 1, klon
	1256	k = ldepar
	1257	kb(i) = ldepar
	1258	ldcum(i) = .TRUE.
	1259	IF (new_deh) THEN
	1260	zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
[524]	1261	ELSE
[1992]	1262	zdeh(i, k) = rcpd(t(i,k-1)-t(i,k)) - rd0.5*(t(i,k-1)+t(i,k))/paprs( &
	1263	i, k)(pplay(i,k-1)-pplay(i,k)) + rlvtt(zqs(i,k-1)-zqs(i,k))
	1264	END IF
	1265	zdeh(i, k) = zdeh(i, k)*0.5
	1266	END DO
	1267	DO k = ldepar + 1, klev
[524]	1268	DO i = 1, klon
[1992]	1269	IF (new_deh) THEN
	1270	zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
	1271	ELSE
	1272	zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
	1273	rd0.5(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
	1274	rlvtt*(zqs(i,k-1)-zqs(i,k))
	1275	END IF
	1276	END DO
	1277	END DO
	1278	END IF
	1279
	1280	! -----Chercher le sommet du nuage
	1281	! -----Calculer la convergence de l'humidite (en kg/m**2 a un facteur
	1282	! -----psolpa/RG pres) du bas jusqu'au sommet du nuage.
	1283	! -----Calculer la convergence virtuelle pour que toute la maille
	1284	! -----deviennt nuageuse (du bas jusqu'au sommet du nuage)
	1285
	1286	DO i = 1, klon
	1287	nuage(i) = .TRUE.
	1288	zconv(i) = 0.0
	1289	zvirt(i) = 0.0
	1290	kh(i) = -999
	1291	END DO
	1292	DO k = 1, klev
	1293	DO i = 1, klon
	1294	IF (k>=kb(i) .AND. ldcum(i)) THEN
	1295	nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
	1296	IF (nuage(i)) THEN
	1297	kh(i) = k
	1298	zconv(i) = zconv(i) + conv_q(i, k)dtime(paprs(i,k)-paprs(i,k+1))
	1299	zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k) &
	1300	-paprs(i,k+1))
	1301	END IF
	1302	END IF
	1303	END DO
	1304	END DO
	1305
	1306	DO i = 1, klon
	1307	todo(i) = ldcum(i) .AND. kh(i) > kb(i) .AND. zconv(i) > 0.0
	1308	END DO
	1309
	1310	kbmin = klev
	1311	kbmax = 0
	1312	khmax = 0
	1313	DO i = 1, klon
	1314	IF (todo(i)) THEN
	1315	kbmin = min(kbmin, kb(i))
	1316	kbmax = max(kbmax, kb(i))
	1317	khmax = max(khmax, kh(i))
	1318	END IF
	1319	END DO
	1320
	1321	! -----Calculer la surface couverte par le nuage
	1322
	1323	DO i = 1, klon
	1324	IF (todo(i)) THEN
	1325	zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
	1326	END IF
	1327	END DO
	1328
	1329	! -----Calculs essentiels:
	1330
	1331	DO i = 1, klon
	1332	IF (todo(i)) THEN
	1333	zcond(i) = 0.0
	1334	END IF
	1335	END DO
	1336	DO k = kbmin, khmax
	1337	DO i = 1, klon
	1338	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1339	zvar = zdeh(i, k)/(1.+zdqs(i,k))
	1340	d_t(i, k) = zvar*zfrac(i)/rcpd
	1341	d_q(i, k) = (zvarzdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))zfrac(i) - &
	1342	conv_q(i, k)*dtime
	1343	zcond(i) = zcond(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
	1344	rneb(i, k) = zfrac(i)
	1345	END IF
	1346	END DO
	1347	END DO
	1348
	1349	DO i = 1, klon
	1350	IF (todo(i) .AND. zcond(i)<0.0) THEN
	1351	PRINT *, 'WARNING: cond. negative (Kuo) ', i, kb(i), kh(i), zcond(i)
	1352	zcond(i) = 0.0
	1353	DO k = kb(i), kh(i)
	1354	d_t(i, k) = 0.0
	1355	d_q(i, k) = 0.0
	1356	END DO
	1357	todo(i) = .FALSE. ! effort totalement perdu
	1358	END IF
	1359	END DO
	1360
	1361	! =====
	1362	! Une fois que la condensation a lieu, on doit construire un
	1363	! "modele nuageux" pour partager la condensation entre l'eau
	1364	! liquide nuageuse et la precipitation (leur rapport toliq
	1365	! est calcule selon l'epaisseur nuageuse). Je suppose que
	1366	! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
	1367	! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
	1368	! lineaire entre dpmin et dpmax).
	1369	! =====
	1370	DO i = 1, klon
	1371	IF (todo(i)) THEN
	1372	toliq(i) = tomax - ((paprs(i,kb(i))-paprs(i,kh(i)+1))/paprs(i,1)-dpmin) &
	1373	*(tomax-tomin)/(dpmax-dpmin)
	1374	toliq(i) = max(tomin, min(tomax,toliq(i)))
	1375	IF (pplay(i,kh(i))/paprs(i,1)<=deep_sig) toliq(i) = deep_to
	1376	IF (old_tau) toliq(i) = 1.0
	1377	END IF
	1378	END DO
	1379	! =====
	1380	! On doit aussi determiner la distribution verticale de
	1381	! l'eau nuageuse. Plusieurs options sont proposees:
	1382
	1383	! (0) La condensation precipite integralement (toliq ne sera
	1384	! pas utilise).
	1385	! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
	1386	! a la vapeur d'eau locale.
	1387	! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
	1388	! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
	1389	! est effectivement diminuee pendant le processus d'ajustement.
	1390	! (4) Elle est en fonction (lineaire ou exponentielle) de la
	1391	! distance (epaisseur en pression) avec le niveau k1 (la couche
	1392	! k1 n'aura donc pas d'eau liquide).
	1393	! =====
	1394
	1395	IF (opt_cld==0) THEN
	1396
	1397	DO i = 1, klon
	1398	IF (todo(i)) zrfl(i) = zcond(i)/dtime
	1399	END DO
	1400
	1401	ELSE IF (opt_cld==1) THEN
	1402
	1403	DO i = 1, klon
	1404	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
	1405	END DO
	1406	DO k = kbmin, khmax
[524]	1407	DO i = 1, klon
[1992]	1408	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1409	zvapo(i) = zvapo(i) + (q(i,k)+d_q(i,k))*(paprs(i,k)-paprs(i,k+1))/ &
	1410	rg
	1411	END IF
	1412	END DO
	1413	END DO
	1414	DO i = 1, klon
[524]	1415	IF (todo(i)) THEN
[1992]	1416	zrapp(i) = toliq(i)*zcond(i)/zvapo(i)
	1417	zrapp(i) = max(0., min(1.,zrapp(i)))
	1418	END IF
	1419	END DO
	1420	DO k = kbmin, khmax
[524]	1421	DO i = 1, klon
[1992]	1422	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1423	d_ql(i, k) = zrapp(i)*(q(i,k)+d_q(i,k))
	1424	END IF
	1425	END DO
	1426	END DO
	1427	DO i = 1, klon
[524]	1428	IF (todo(i)) THEN
[1992]	1429	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	1430	END IF
	1431	END DO
	1432
	1433	ELSE IF (opt_cld==2) THEN
	1434
	1435	DO i = 1, klon
	1436	IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
	1437	END DO
	1438	DO k = kbmin, khmax
[524]	1439	DO i = 1, klon
[1992]	1440	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1441	zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg
	1442	END IF
	1443	END DO
	1444	END DO
	1445	DO k = kbmin, khmax
	1446	DO i = 1, klon
	1447	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1448	d_ql(i, k) = toliq(i)*zcond(i)/zvapo(i)
	1449	END IF
	1450	END DO
	1451	END DO
	1452	DO i = 1, klon
[524]	1453	IF (todo(i)) THEN
[1992]	1454	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	1455	END IF
	1456	END DO
	1457
	1458	ELSE IF (opt_cld==3) THEN
	1459
	1460	DO i = 1, klon
	1461	IF (todo(i)) THEN
	1462	zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
	1463	END IF
	1464	END DO
	1465	DO k = kbmin, khmax
[524]	1466	DO i = 1, klon
[1992]	1467	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
	1468	zvapo(i) = zvapo(i) + max(0.0, -d_q(i,k))*(paprs(i,k)-paprs(i,k+1)) &
	1469	/rg
	1470	END IF
	1471	END DO
	1472	END DO
	1473	DO k = kbmin, khmax
[524]	1474	DO i = 1, klon
[1992]	1475	IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i) .AND. zvapo(i)>0.0) THEN
	1476	d_ql(i, k) = d_ql(i, k) + toliq(i)zcond(i)/zvapo(i)max(0.0, -d_q( &
	1477	i,k))
	1478	END IF
	1479	END DO
	1480	END DO
	1481	DO i = 1, klon
[524]	1482	IF (todo(i)) THEN
[1992]	1483	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	1484	END IF
	1485	END DO
	1486
	1487	ELSE IF (opt_cld==4) THEN
	1488
	1489	nexpo = 3
	1490	! cc nexpo = 1 ! distribution lineaire
	1491
	1492	DO i = 1, klon
	1493	IF (todo(i)) THEN
	1494	zvapo(i) = 0.0 ! quantite integrale de masse (avec ponderation)
	1495	END IF
	1496	END DO
	1497	DO k = kbmin, khmax
[524]	1498	DO i = 1, klon
[1992]	1499	IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
	1500	zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg*(pplay(i,kb(i))- &
	1501	pplay(i,k))**nexpo
	1502	END IF
	1503	END DO
	1504	END DO
	1505	DO k = kbmin, khmax
[524]	1506	DO i = 1, klon
[1992]	1507	IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
	1508	d_ql(i, k) = d_ql(i, k) + toliq(i)zcond(i)/zvapo(i)(pplay(i,kb(i) &
	1509	)-pplay(i,k))**nexpo
	1510	END IF
	1511	END DO
	1512	END DO
	1513	DO i = 1, klon
[524]	1514	IF (todo(i)) THEN
[1992]	1515	zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
	1516	END IF
	1517	END DO
	1518
	1519	ELSE ! valeur non-prevue pour opt_cld
	1520
	1521	PRINT *, 'opt_cld est faux:', opt_cld
	1522	CALL abort
	1523
	1524	END IF ! fin de opt_cld
	1525
	1526	! L'eau precipitante peut etre re-evaporee:
	1527
	1528	IF (evap_prec .AND. kbmax>=2) THEN
	1529	DO k = kbmax, 1, -1
	1530	DO i = 1, klon
	1531	IF (todo(i) .AND. k<=(kb(i)-1) .AND. zrfl(i)>0.0) THEN
	1532	zqev = max(0.0, (zqs(i,k)-q(i,k))*zfrac(i))
	1533	zqevt = coef_eva(1.0-q(i,k)/zqs(i,k))sqrt(zrfl(i))* &
	1534	(paprs(i,k)-paprs(i,k+1))/pplay(i, k)t(i, k)rd/rg
	1535	zqevt = max(0.0, min(zqevt,zrfl(i)))rgdtime/ &
	1536	(paprs(i,k)-paprs(i,k+1))
	1537	zqev = min(zqev, zqevt)
	1538	zrfln = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime
	1539	d_q(i, k) = -(zrfln-zrfl(i))(rg/(paprs(i,k)-paprs(i,k+1)))dtime
	1540	d_t(i, k) = (zrfln-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
	1541	k+1)))dtimerlvtt/rcpd
	1542	zrfl(i) = zrfln
	1543	END IF
	1544	END DO
	1545	END DO
	1546	END IF
	1547
	1548	! La temperature de la premiere couche determine la pluie ou la neige:
	1549
	1550	DO i = 1, klon
	1551	IF (todo(i)) THEN
	1552	IF (t(i,1)>rtt) THEN
	1553	rain(i) = rain(i) + zrfl(i)
[524]	1554	ELSE
[1992]	1555	snow(i) = snow(i) + zrfl(i)
	1556	END IF
	1557	END IF
	1558	END DO
	1559
	1560	RETURN
	1561	END SUBROUTINE conkuo
	1562	SUBROUTINE kuofcl(pt, pq, pg, pp, ldcum, kcbot)
	1563	USE dimphy
	1564	IMPLICIT NONE
	1565	! ======================================================================
	1566	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
	1567	! adaptation du code de Tiedtke du ECMWF
	1568	! Objet: calculer le niveau de convection libre
	1569	! (FCL: Free Convection Level)
	1570	! ======================================================================
	1571	! Arguments:
	1572	! pt---input-R- temperature (K)
	1573	! pq---input-R- vapeur d'eau (kg/kg)
	1574	! pg---input-R- geopotentiel (g*z ou z est en metre)
	1575	! pp---input-R- pression (Pa)
	1576
	1577	! LDCUM---output-L- Y-t-il la convection
	1578	! kcbot---output-I- Niveau du bas de la convection
	1579	! ======================================================================
	1580	! ym#include "dimensions.h"
	1581	! ym#include "dimphy.h"
	1582	include "YOMCST.h"
	1583	include "YOETHF.h"
	1584
	1585	REAL pt(klon, klev), pq(klon, klev), pg(klon, klev), pp(klon, klev)
	1586	INTEGER kcbot(klon)
	1587	LOGICAL ldcum(klon)
	1588
	1589	REAL ztu(klon, klev), zqu(klon, klev), zlu(klon, klev)
	1590	REAL zqold(klon), zbuo
	1591	INTEGER is, i, k
	1592
	1593	! klab=1: on est sous le nuage convectif
	1594	! klab=2: le bas du nuage convectif
	1595	! klab=0: autres couches
	1596	INTEGER klab(klon, klev)
	1597
	1598	! quand lflag=.true., on est sous le nuage, il faut donc appliquer
	1599	! le processus d'elevation.
	1600	LOGICAL lflag(klon)
	1601
	1602	DO k = 1, klev
	1603	DO i = 1, klon
	1604	ztu(i, k) = pt(i, k)
	1605	zqu(i, k) = pq(i, k)
	1606	zlu(i, k) = 0.0
	1607	klab(i, k) = 0
	1608	END DO
	1609	END DO
	1610	! ----------------------------------------------------------------------
	1611	DO i = 1, klon
	1612	klab(i, 1) = 1
	1613	kcbot(i) = 2
	1614	ldcum(i) = .FALSE.
	1615	END DO
	1616
	1617	DO k = 2, klev - 1
	1618
	1619	is = 0
	1620	DO i = 1, klon
	1621	IF (klab(i,k-1)==1) is = is + 1
	1622	lflag(i) = .FALSE.
	1623	IF (klab(i,k-1)==1) lflag(i) = .TRUE.
	1624	END DO
	1625	IF (is==0) GO TO 290
	1626
	1627	! on eleve le parcel d'air selon l'adiabatique sec
	1628
	1629	DO i = 1, klon
[524]	1630	IF (lflag(i)) THEN
[1992]	1631	zqu(i, k) = zqu(i, k-1)
	1632	ztu(i, k) = ztu(i, k-1) + (pg(i,k-1)-pg(i,k))/rcpd
	1633	zbuo = ztu(i, k)(1.+retvzqu(i,k)) - pt(i, k)(1.+retvpq(i,k)) + &
	1634	0.5
	1635	IF (zbuo>0.) klab(i, k) = 1
	1636	zqold(i) = zqu(i, k)
[524]	1637	END IF
[1992]	1638	END DO
[524]	1639
[1992]	1640	! on calcule la condensation eventuelle
	1641
	1642	CALL adjtq(pp(1,k), ztu(1,k), zqu(1,k), lflag, 1)
	1643
	1644	! s'il y a la condensation et la "buoyancy" force est positive
	1645	! c'est bien le bas de la tour de convection
	1646
	1647	DO i = 1, klon
	1648	IF (lflag(i) .AND. zqu(i,k)/=zqold(i)) THEN
	1649	klab(i, k) = 2
	1650	zlu(i, k) = zlu(i, k) + zqold(i) - zqu(i, k)
	1651	zbuo = ztu(i, k)(1.+retvzqu(i,k)) - pt(i, k)(1.+retvpq(i,k)) + &
	1652	0.5
	1653	IF (zbuo>0.) THEN
	1654	kcbot(i) = k
	1655	ldcum(i) = .TRUE.
	1656	END IF
	1657	END IF
	1658	END DO
	1659
	1660	290 END DO
	1661
	1662	RETURN
	1663	END SUBROUTINE kuofcl
	1664	SUBROUTINE adjtq(pp, pt, pq, ldflag, kcall)
	1665	USE dimphy
	1666	IMPLICIT NONE
	1667	! ======================================================================
	1668	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
	1669	! adaptation du code de Tiedtke du ECMWF
	1670	! Objet: ajustement entre T et Q
	1671	! ======================================================================
	1672	! Arguments:
	1673	! pp---input-R- pression (Pa)
	1674	! pt---input/output-R- temperature (K)
	1675	! pq---input/output-R- vapeur d'eau (kg/kg)
	1676	! ======================================================================
	1677	! TO PRODUCE T,Q AND L VALUES FOR CLOUD ASCENT
	1678
	1679	! NOTE: INPUT PARAMETER KCALL DEFINES CALCULATION AS
	1680	! KCALL=0 ENV. T AND QS INCUINI
	1681	! KCALL=1 CONDENSATION IN UPDRAFTS (E.G. CUBASE, CUASC)
	1682	! KCALL=2 EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)
	1683
	1684	! ym#include "dimensions.h"
	1685	! ym#include "dimphy.h"
	1686	include "YOMCST.h"
	1687
	1688	REAL pt(klon), pq(klon), pp(klon)
	1689	LOGICAL ldflag(klon)
	1690	INTEGER kcall
	1691
	1692	REAL t_coup
	1693	PARAMETER (t_coup=234.0)
	1694
	1695	REAL zcond(klon), zcond1
	1696	REAL zdelta, zcvm5, zldcp, zqsat, zcor, zdqsat
	1697	INTEGER is, i
	1698	include "YOETHF.h"
	1699	include "FCTTRE.h"
	1700
	1701	DO i = 1, klon
	1702	zcond(i) = 0.0
	1703	END DO
	1704
	1705	DO i = 1, klon
	1706	IF (ldflag(i)) THEN
	1707	zdelta = max(0., sign(1.,rtt-pt(i)))
	1708	zldcp = rlvtt(1.-zdelta) + zdeltarlstt
	1709	zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
	1710	IF (thermcep) THEN
	1711	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1712	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
	1713	zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
	1714	zqsat = min(0.5, zqsat)
	1715	zcor = 1./(1.-retv*zqsat)
	1716	zqsat = zqsat*zcor
	1717	zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
[524]	1718	ELSE
[1992]	1719	IF (pt(i)<t_coup) THEN
	1720	zqsat = qsats(pt(i))/pp(i)
	1721	zdqsat = dqsats(pt(i), zqsat)
	1722	ELSE
	1723	zqsat = qsatl(pt(i))/pp(i)
	1724	zdqsat = dqsatl(pt(i), zqsat)
	1725	END IF
	1726	END IF
	1727	zcond(i) = (pq(i)-zqsat)/(1.+zdqsat)
	1728	IF (kcall==1) zcond(i) = max(zcond(i), 0.)
	1729	IF (kcall==2) zcond(i) = min(zcond(i), 0.)
	1730	pt(i) = pt(i) + zldcp*zcond(i)
	1731	pq(i) = pq(i) - zcond(i)
	1732	END IF
	1733	END DO
	1734
	1735	is = 0
	1736	DO i = 1, klon
	1737	IF (zcond(i)/=0.) is = is + 1
	1738	END DO
	1739	IF (is==0) GO TO 230
	1740
	1741	DO i = 1, klon
	1742	IF (ldflag(i) .AND. zcond(i)/=0.) THEN
	1743	zdelta = max(0., sign(1.,rtt-pt(i)))
	1744	zldcp = rlvtt(1.-zdelta) + zdeltarlstt
	1745	zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
	1746	IF (thermcep) THEN
	1747	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1748	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
	1749	zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
	1750	zqsat = min(0.5, zqsat)
	1751	zcor = 1./(1.-retv*zqsat)
	1752	zqsat = zqsat*zcor
	1753	zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
[524]	1754	ELSE
[1992]	1755	IF (pt(i)<t_coup) THEN
	1756	zqsat = qsats(pt(i))/pp(i)
	1757	zdqsat = dqsats(pt(i), zqsat)
	1758	ELSE
	1759	zqsat = qsatl(pt(i))/pp(i)
	1760	zdqsat = dqsatl(pt(i), zqsat)
	1761	END IF
	1762	END IF
	1763	zcond1 = (pq(i)-zqsat)/(1.+zdqsat)
	1764	pt(i) = pt(i) + zldcp*zcond1
	1765	pq(i) = pq(i) - zcond1
	1766	END IF
	1767	END DO
	1768
	1769	230 CONTINUE
	1770	RETURN
	1771	END SUBROUTINE adjtq
	1772	SUBROUTINE fiajh(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, &
	1773	ibas, itop)
	1774	USE dimphy
	1775	IMPLICIT NONE
	1776
	1777	! Ajustement humide (Schema de convection de Manabe)
	1778	! .
	1779	! ym#include "dimensions.h"
	1780	! ym#include "dimphy.h"
	1781	include "YOMCST.h"
	1782
	1783	! Arguments:
	1784
	1785	REAL dtime ! intervalle du temps (s)
	1786	REAL t(klon, klev) ! temperature (K)
	1787	REAL q(klon, klev) ! humidite specifique (kg/kg)
	1788	REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
	1789	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	1790
	1791	REAL d_t(klon, klev) ! incrementation pour la temperature
	1792	REAL d_q(klon, klev) ! incrementation pour vapeur d'eau
	1793	REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
	1794	REAL rneb(klon, klev) ! fraction nuageuse
	1795
	1796	REAL rain(klon) ! variable non utilisee
	1797	REAL snow(klon) ! variable non utilisee
	1798	INTEGER ibas(klon) ! variable non utilisee
	1799	INTEGER itop(klon) ! variable non utilisee
	1800
	1801	REAL t_coup
	1802	PARAMETER (t_coup=234.0)
	1803	REAL seuil_vap
	1804	PARAMETER (seuil_vap=1.0E-10)
	1805
	1806	! Variables locales:
	1807
	1808	INTEGER i, k
	1809	INTEGER k1, k1p, k2, k2p
	1810	LOGICAL itest(klon)
	1811	REAL delta_q(klon, klev)
	1812	REAL cp_new_t(klev)
	1813	REAL cp_delta_t(klev)
	1814	REAL new_qb(klev)
	1815	REAL v_cptj(klev), v_cptjk1, v_ssig
	1816	REAL v_cptt(klon, klev), v_p, v_t
	1817	REAL v_qs(klon, klev), v_qsd(klon, klev)
	1818	REAL zq1(klon), zq2(klon)
	1819	REAL gamcpdz(klon, 2:klev)
	1820	REAL zdp, zdpm
	1821
	1822	REAL zsat ! sur-saturation
	1823	REAL zflo ! flotabilite
	1824
	1825	REAL local_q(klon, klev), local_t(klon, klev)
	1826
	1827	REAL zdelta, zcor, zcvm5
	1828
	1829	include "YOETHF.h"
	1830	include "FCTTRE.h"
	1831
	1832	DO k = 1, klev
	1833	DO i = 1, klon
	1834	local_q(i, k) = q(i, k)
	1835	local_t(i, k) = t(i, k)
	1836	rneb(i, k) = 0.0
	1837	d_ql(i, k) = 0.0
	1838	d_t(i, k) = 0.0
	1839	d_q(i, k) = 0.0
	1840	END DO
	1841	END DO
	1842	DO i = 1, klon
	1843	rain(i) = 0.0
	1844	snow(i) = 0.0
	1845	ibas(i) = 0
	1846	itop(i) = 0
	1847	END DO
	1848
	1849	! Calculer v_qs et v_qsd:
	1850
	1851	DO k = 1, klev
	1852	DO i = 1, klon
	1853	v_cptt(i, k) = rcpd*local_t(i, k)
	1854	v_t = local_t(i, k)
	1855	v_p = pplay(i, k)
	1856
	1857	IF (thermcep) THEN
	1858	zdelta = max(0., sign(1.,rtt-v_t))
	1859	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1860	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
	1861	v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
	1862	v_qs(i, k) = min(0.5, v_qs(i,k))
	1863	zcor = 1./(1.-retv*v_qs(i,k))
	1864	v_qs(i, k) = v_qs(i, k)*zcor
	1865	v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
	1866	ELSE
	1867	IF (v_t<t_coup) THEN
	1868	v_qs(i, k) = qsats(v_t)/v_p
	1869	v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
	1870	ELSE
	1871	v_qs(i, k) = qsatl(v_t)/v_p
	1872	v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
	1873	END IF
	1874	END IF
	1875	END DO
	1876	END DO
	1877
	1878	! Calculer Gamma * Cp * dz: (gamm est le gradient critique)
	1879
	1880	DO k = 2, klev
	1881	DO i = 1, klon
	1882	zdp = paprs(i, k) - paprs(i, k+1)
	1883	zdpm = paprs(i, k-1) - paprs(i, k)
	1884	gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)(v_cptt(i,k-1)zdpm+ &
	1885	v_cptt(i,k)zdp)+rlvtt/(zdpm+zdp)(v_qs(i,k-1)*zdpm+ &
	1886	v_qs(i,k)zdp))(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
	1887	k-1)zdpm+v_qsd(i,k)zdp)/(zdpm+zdp))
	1888	END DO
	1889	END DO
	1890
	1891	! ------------------------------------ modification des profils instables
	1892	DO i = 1, klon
	1893	itest(i) = .FALSE.
	1894
	1895	k1 = 0
	1896	k2 = 1
	1897
	1898	810 CONTINUE ! chercher k1, le bas de la colonne
	1899	k2 = k2 + 1
	1900	IF (k2>klev) GO TO 9999
	1901	zflo = v_cptt(i, k2-1) - v_cptt(i, k2) - gamcpdz(i, k2)
	1902	zsat = (local_q(i,k2-1)-v_qs(i,k2-1))*(paprs(i,k2-1)-paprs(i,k2)) + &
	1903	(local_q(i,k2)-v_qs(i,k2))*(paprs(i,k2)-paprs(i,k2+1))
	1904	IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 810
	1905	k1 = k2 - 1
	1906	itest(i) = .TRUE.
	1907
	1908	820 CONTINUE ! chercher k2, le haut de la colonne
	1909	IF (k2==klev) GO TO 821
	1910	k2p = k2 + 1
	1911	zsat = zsat + (paprs(i,k2p)-paprs(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
	1912	zflo = v_cptt(i, k2p-1) - v_cptt(i, k2p) - gamcpdz(i, k2p)
	1913	IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 821
	1914	k2 = k2p
	1915	GO TO 820
	1916	821 CONTINUE
	1917
	1918	! ------------------------------------------------------ ajustement local
	1919	830 CONTINUE ! ajustement proprement dit
	1920	v_cptj(k1) = 0.0
	1921	zdp = paprs(i, k1) - paprs(i, k1+1)
	1922	v_cptjk1 = ((1.0+v_qsd(i,k1))(v_cptt(i,k1)+v_cptj(k1))+rlvtt(local_q(i, &
	1923	k1)-v_qs(i,k1)))*zdp
	1924	v_ssig = zdp*(1.0+v_qsd(i,k1))
	1925
	1926	k1p = k1 + 1
	1927	DO k = k1p, k2
	1928	zdp = paprs(i, k) - paprs(i, k+1)
	1929	v_cptj(k) = v_cptj(k-1) + gamcpdz(i, k)
	1930	v_cptjk1 = v_cptjk1 + zdp((1.0+v_qsd(i,k))(v_cptt(i, &
	1931	k)+v_cptj(k))+rlvtt*(local_q(i,k)-v_qs(i,k)))
	1932	v_ssig = v_ssig + zdp*(1.0+v_qsd(i,k))
	1933	END DO
	1934
	1935	DO k = k1, k2
	1936	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
	1937	cp_delta_t(k) = cp_new_t(k) - v_cptt(i, k)
	1938	new_qb(k) = v_qs(i, k) + v_qsd(i, k)*cp_delta_t(k)/rlvtt
	1939	local_q(i, k) = new_qb(k)
	1940	local_t(i, k) = cp_new_t(k)/rcpd
	1941	END DO
	1942
	1943	! --------------------------------------------------- sondage vers le bas
	1944	! -- on redefinit les variables prognostiques dans
	1945	! -- la colonne qui vient d'etre ajustee
	1946
	1947	DO k = k1, k2
	1948	v_cptt(i, k) = rcpd*local_t(i, k)
	1949	v_t = local_t(i, k)
	1950	v_p = pplay(i, k)
	1951
	1952	IF (thermcep) THEN
	1953	zdelta = max(0., sign(1.,rtt-v_t))
	1954	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	1955	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
	1956	v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
	1957	v_qs(i, k) = min(0.5, v_qs(i,k))
	1958	zcor = 1./(1.-retv*v_qs(i,k))
	1959	v_qs(i, k) = v_qs(i, k)*zcor
	1960	v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
	1961	ELSE
	1962	IF (v_t<t_coup) THEN
	1963	v_qs(i, k) = qsats(v_t)/v_p
	1964	v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
	1965	ELSE
	1966	v_qs(i, k) = qsatl(v_t)/v_p
	1967	v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
	1968	END IF
	1969	END IF
	1970	END DO
	1971	DO k = 2, klev
	1972	zdpm = paprs(i, k-1) - paprs(i, k)
	1973	zdp = paprs(i, k) - paprs(i, k+1)
	1974	gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)(v_cptt(i,k-1)zdpm+ &
	1975	v_cptt(i,k)zdp)+rlvtt/(zdpm+zdp)(v_qs(i,k-1)*zdpm+ &
	1976	v_qs(i,k)zdp))(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
	1977	k-1)zdpm+v_qsd(i,k)zdp)/(zdpm+zdp))
	1978	END DO
	1979
	1980	! Verifier si l'on peut etendre la colonne vers le bas
	1981
	1982	IF (k1==1) GO TO 841 ! extension echouee
	1983	zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
	1984	zsat = (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1)) + &
	1985	(local_q(i,k1)-v_qs(i,k1))*(paprs(i,k1)-paprs(i,k1+1))
	1986	IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 841 ! extension echouee
	1987
	1988	840 CONTINUE
	1989	k1 = k1 - 1
	1990	IF (k1==1) GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	1991	zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1))
	1992	zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
	1993	IF (zflo>0.0 .AND. zsat>0.0) THEN
	1994	GO TO 840
	1995	ELSE
	1996	GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	1997	END IF
	1998	841 CONTINUE
	1999
	2000	GO TO 810 ! chercher d'autres blocks en haut
	2001
	2002	9999 END DO ! boucle sur tous les points
	2003	! -----------------------------------------------------------------------
	2004
	2005	! Determiner la fraction nuageuse (hypothese: la nebulosite a lieu
	2006	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
	2007
	2008	DO k = 1, klev
	2009	DO i = 1, klon
	2010	IF (itest(i)) THEN
	2011	delta_q(i, k) = local_q(i, k) - q(i, k)
	2012	IF (delta_q(i,k)<0.) rneb(i, k) = 1.0
	2013	END IF
	2014	END DO
	2015	END DO
	2016
	2017	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
	2018	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
	2019	! diminue et d'une maniere proportionnelle a cet diminution):
	2020
	2021	DO i = 1, klon
	2022	IF (itest(i)) THEN
	2023	zq1(i) = 0.0
	2024	zq2(i) = 0.0
	2025	END IF
	2026	END DO
	2027	DO k = 1, klev
	2028	DO i = 1, klon
	2029	IF (itest(i)) THEN
	2030	zdp = paprs(i, k) - paprs(i, k+1)
	2031	zq1(i) = zq1(i) - delta_q(i, k)*zdp
	2032	zq2(i) = zq2(i) - min(0.0, delta_q(i,k))*zdp
	2033	END IF
	2034	END DO
	2035	END DO
	2036	DO k = 1, klev
	2037	DO i = 1, klon
	2038	IF (itest(i)) THEN
	2039	IF (zq2(i)/=0.0) d_ql(i, k) = -min(0.0, delta_q(i,k))*zq1(i)/zq2(i)
	2040	END IF
	2041	END DO
	2042	END DO
	2043
	2044	DO k = 1, klev
	2045	DO i = 1, klon
	2046	local_q(i, k) = max(local_q(i,k), seuil_vap)
	2047	END DO
	2048	END DO
	2049
	2050	DO k = 1, klev
	2051	DO i = 1, klon
	2052	d_t(i, k) = local_t(i, k) - t(i, k)
	2053	d_q(i, k) = local_q(i, k) - q(i, k)
	2054	END DO
	2055	END DO
	2056
	2057	RETURN
	2058	END SUBROUTINE fiajh
	2059	SUBROUTINE fiajc(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
	2060	rain, snow, ibas, itop)
	2061	USE dimphy
	2062	IMPLICIT NONE
	2063
	2064	! ym#include "dimensions.h"
	2065	! ym#include "dimphy.h"
	2066	include "YOMCST.h"
	2067
	2068	! Options:
	2069
	2070	INTEGER plb ! niveau de depart pour la convection
	2071	PARAMETER (plb=4)
	2072
	2073	! Mystere: cette option n'est pas innocente pour les resultats !
	2074	! Qui peut resoudre ce mystere ? (Z.X.Li mars 1995)
	2075	LOGICAL vector ! calcul vectorise
	2076	PARAMETER (vector=.FALSE.)
	2077
	2078	REAL t_coup
	2079	PARAMETER (t_coup=234.0)
	2080
	2081	! Arguments:
	2082
	2083	REAL q(klon, klev) ! humidite specifique (kg/kg)
	2084	REAL t(klon, klev) ! temperature (K)
	2085	REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
	2086	REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
	2087	REAL dtime ! intervalle du temps (s)
	2088	REAL conv_q(klon, klev) ! taux de convergence de l'humidite
	2089	REAL rneb(klon, klev) ! fraction nuageuse
	2090	REAL d_q(klon, klev) ! incrementaion pour la vapeur d'eau
	2091	REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
	2092	REAL d_t(klon, klev) ! incrementation pour la temperature
	2093	REAL rain(klon) ! variable non-utilisee
	2094	REAL snow(klon) ! variable non-utilisee
	2095	INTEGER itop(klon) ! variable non-utilisee
	2096	INTEGER ibas(klon) ! variable non-utilisee
	2097
	2098	INTEGER kh(klon), i, k
	2099	LOGICAL nuage(klon), test(klon, klev)
	2100	REAL zconv(klon), zdeh(klon, klev), zvirt(klon)
	2101	REAL zdqs(klon, klev), zqs(klon, klev)
	2102	REAL ztt, zvar, zfrac(klon)
	2103	REAL zq1(klon), zq2(klon)
	2104	REAL zdelta, zcor, zcvm5
	2105
	2106	include "YOETHF.h"
	2107	include "FCTTRE.h"
	2108
	2109	! Initialiser les sorties:
	2110
	2111	DO k = 1, klev
	2112	DO i = 1, klon
	2113	rneb(i, k) = 0.0
	2114	d_ql(i, k) = 0.0
	2115	d_t(i, k) = 0.0
	2116	d_q(i, k) = 0.0
	2117	END DO
	2118	END DO
	2119	DO i = 1, klon
	2120	itop(i) = 0
	2121	ibas(i) = 0
	2122	rain(i) = 0.0
	2123	snow(i) = 0.0
	2124	END DO
	2125
	2126	! Calculer Qs et L/Cp * dQs/dT:
	2127
	2128	DO k = 1, klev
	2129	DO i = 1, klon
	2130	ztt = t(i, k)
	2131	IF (thermcep) THEN
	2132	zdelta = max(0., sign(1.,rtt-ztt))
	2133	zcvm5 = r5lesrlvtt(1.-zdelta) + zdeltar5iesrlstt
	2134	zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
	2135	zqs(i, k) = r2es*foeew(ztt, zdelta)/pplay(i, k)
	2136	zqs(i, k) = min(0.5, zqs(i,k))
	2137	zcor = 1./(1.-retv*zqs(i,k))
	2138	zqs(i, k) = zqs(i, k)*zcor
	2139	zdqs(i, k) = foede(ztt, zdelta, zcvm5, zqs(i,k), zcor)
	2140	ELSE
	2141	IF (ztt<t_coup) THEN
	2142	zqs(i, k) = qsats(ztt)/pplay(i, k)
	2143	zdqs(i, k) = dqsats(ztt, zqs(i,k))
	2144	ELSE
	2145	zqs(i, k) = qsatl(ztt)/pplay(i, k)
	2146	zdqs(i, k) = dqsatl(ztt, zqs(i,k))
	2147	END IF
	2148	END IF
	2149	END DO
	2150	END DO
	2151
	2152	! Determiner la difference de l'energie totale saturee:
	2153
	2154	DO i = 1, klon
	2155	k = plb
	2156	zdeh(i, k) = rcpd(t(i,k-1)-t(i,k)) - rd0.5*(t(i,k-1)+t(i,k))/paprs(i, k &
	2157	)(pplay(i,k-1)-pplay(i,k)) + rlvtt(zqs(i,k-1)-zqs(i,k))
	2158	zdeh(i, k) = zdeh(i, k)*0.5 ! on prend la moitie
	2159	END DO
	2160	DO k = plb + 1, klev
	2161	DO i = 1, klon
	2162	zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
	2163	rd0.5(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
	2164	rlvtt*(zqs(i,k-1)-zqs(i,k))
	2165	END DO
	2166	END DO
	2167
	2168	! Determiner le sommet du nuage selon l'instabilite
	2169	! Calculer les convergences d'humidite (reelle et virtuelle)
	2170
	2171	DO i = 1, klon
	2172	nuage(i) = .TRUE.
	2173	zconv(i) = 0.0
	2174	zvirt(i) = 0.0
	2175	kh(i) = -999
	2176	END DO
	2177	DO k = plb, klev
	2178	DO i = 1, klon
	2179	nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
	2180	IF (nuage(i)) THEN
	2181	kh(i) = k
	2182	zconv(i) = zconv(i) + conv_q(i, k)dtime(paprs(i,k)-paprs(i,k+1))
	2183	zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k)- &
	2184	paprs(i,k+1))
	2185	END IF
	2186	END DO
	2187	END DO
	2188
	2189	IF (vector) THEN
	2190
	2191
	2192	DO k = plb, klev
[524]	2193	DO i = 1, klon
[1992]	2194	IF (k<=kh(i) .AND. kh(i)>plb .AND. zconv(i)>0.0) THEN
	2195	test(i, k) = .TRUE.
	2196	zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
	2197	ELSE
	2198	test(i, k) = .FALSE.
	2199	END IF
	2200	END DO
	2201	END DO
	2202
	2203	DO k = plb, klev
[524]	2204	DO i = 1, klon
[1992]	2205	IF (test(i,k)) THEN
	2206	zvar = zdeh(i, k)/(1.0+zdqs(i,k))
	2207	d_q(i, k) = (zvarzdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))zfrac(i) - &
	2208	conv_q(i, k)*dtime
	2209	d_t(i, k) = zvar*zfrac(i)/rcpd
	2210	END IF
	2211	END DO
	2212	END DO
	2213
	2214	DO i = 1, klon
	2215	zq1(i) = 0.0
	2216	zq2(i) = 0.0
	2217	END DO
	2218	DO k = plb, klev
[524]	2219	DO i = 1, klon
[1992]	2220	IF (test(i,k)) THEN
	2221	IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
	2222	zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
	2223	zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
	2224	END IF
	2225	END DO
	2226	END DO
	2227
	2228	DO k = plb, klev
[524]	2229	DO i = 1, klon
[1992]	2230	IF (test(i,k)) THEN
	2231	IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
	2232	END IF
	2233	END DO
	2234	END DO
	2235
	2236	ELSE ! (.NOT. vector)
	2237
	2238	DO i = 1, klon
	2239	IF (kh(i)>plb .AND. zconv(i)>0.0) THEN
	2240	! cc IF (kh(i).LE.plb) GOTO 999 ! il n'y a pas d'instabilite
	2241	! cc IF (zconv(i).LE.0.0) GOTO 999 ! convergence insuffisante
	2242	zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
	2243	DO k = plb, kh(i)
	2244	zvar = zdeh(i, k)/(1.0+zdqs(i,k))
	2245	d_q(i, k) = (zvarzdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))zfrac(i) - &
	2246	conv_q(i, k)*dtime
	2247	d_t(i, k) = zvar*zfrac(i)/rcpd
	2248	END DO
	2249
	2250	zq1(i) = 0.0
	2251	zq2(i) = 0.0
	2252	DO k = plb, kh(i)
	2253	IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
	2254	zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
	2255	zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
	2256	END DO
	2257	DO k = plb, kh(i)
	2258	IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
	2259	END DO
	2260	END IF
	2261	END DO
	2262
	2263	END IF ! fin de teste sur vector
	2264
	2265	RETURN
	2266	END SUBROUTINE fiajc

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

Download in other formats:

Original Format