Context Navigation

thermcell.F @ 2394

Last change on this file since 2394 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 39.2 KB

Rev	Line
[1299]	1	!
	2	! $Id: thermcell.F 1299 2010-01-20 14:27:21Z acozic $
	3	!
[878]	4	SUBROUTINE calcul_sec(ngrid,nlay,ptimestep
	5	s ,pplay,pplev,pphi,zlev
[542]	6	s ,pu,pv,pt,po
[878]	7	s ,zmax,wmax,zw2,lmix
[542]	8	c s ,pu_therm,pv_therm
	9	s ,r_aspect,l_mix,w2di,tho)
[878]	10
[940]	11	USE dimphy
[542]	12	IMPLICIT NONE
	13
	14	c=======================================================================
	15	c
	16	c Calcul du transport verticale dans la couche limite en presence
	17	c de "thermiques" explicitement representes
	18	c
[878]	19	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
[542]	20	c
[878]	21	c le thermique est supposé homogène et dissipé par mélange avec
	22	c son environnement. la longueur l_mix contrôle l'efficacité du
	23	c mélange
[542]	24	c
[878]	25	c Le calcul du transport des différentes espèces se fait en prenant
[542]	26	c en compte:
	27	c 1. un flux de masse montant
	28	c 2. un flux de masse descendant
	29	c 3. un entrainement
	30	c 4. un detrainement
	31	c
	32	c=======================================================================
	33
	34	c-----------------------------------------------------------------------
	35	c declarations:
	36	c -------------
	37
[878]	38	#include "dimensions.h"
[940]	39	cccc#include "dimphy.h"
[542]	40	#include "YOMCST.h"
	41
	42	c arguments:
	43	c ----------
	44
	45	INTEGER ngrid,nlay,w2di,tho
	46	real ptimestep,l_mix,r_aspect
	47	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
	48	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
	49	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
	50	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
	51	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	52	real pphi(ngrid,nlay)
	53
	54	integer idetr
	55	save idetr
	56	data idetr/3/
[987]	57	c$OMP THREADPRIVATE(idetr)
[542]	58	c local:
	59	c ------
	60
	61	INTEGER ig,k,l,lmaxa(klon),lmix(klon)
	62	real zsortie1d(klon)
	63	c CR: on remplace lmax(klon,klev+1)
	64	INTEGER lmax(klon),lmin(klon),lentr(klon)
	65	real linter(klon)
	66	real zmix(klon), fracazmix(klon)
	67	c RC
[940]	68	real zmax(klon),zw,zw2(klon,klev+1),ztva(klon,klev)
[542]	69
	70	real zlev(klon,klev+1),zlay(klon,klev)
	71	REAL zh(klon,klev),zdhadj(klon,klev)
	72	REAL ztv(klon,klev)
	73	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
	74	REAL wh(klon,klev+1)
	75	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
	76	real zla(klon,klev+1)
	77	real zwa(klon,klev+1)
	78	real zld(klon,klev+1)
[940]	79	! real zwd(klon,klev+1)
[542]	80	real zsortie(klon,klev)
	81	real zva(klon,klev)
	82	real zua(klon,klev)
	83	real zoa(klon,klev)
	84
	85	real zha(klon,klev)
	86	real wa_moy(klon,klev+1)
	87	real fraca(klon,klev+1)
	88	real fracc(klon,klev+1)
	89	real zf,zf2
[878]	90	real thetath2(klon,klev),wth2(klon,klev)
[940]	91	! common/comtherm/thetath2,wth2
[542]	92
	93	real count_time
[940]	94	! integer isplit,nsplit
[542]	95	integer isplit,nsplit,ialt
	96	parameter (nsplit=10)
	97	data isplit/0/
	98	save isplit
[987]	99	c$OMP THREADPRIVATE(isplit)
[878]	100
[542]	101	logical sorties
	102	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
	103	real zpspsk(klon,klev)
	104
	105	c real wmax(klon,klev),wmaxa(klon)
	106	real wmax(klon),wmaxa(klon)
	107	real wa(klon,klev,klev+1)
	108	real wd(klon,klev+1)
	109	real larg_part(klon,klev,klev+1)
	110	real fracd(klon,klev+1)
	111	real xxx(klon,klev+1)
	112	real larg_cons(klon,klev+1)
	113	real larg_detr(klon,klev+1)
	114	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
	115	real pu_therm(klon,klev),pv_therm(klon,klev)
	116	real fm(klon,klev+1),entr(klon,klev)
	117	real fmc(klon,klev+1)
[878]	118
[542]	119	cCR:nouvelles variables
	120	real f_star(klon,klev+1),entr_star(klon,klev)
	121	real entr_star_tot(klon),entr_star2(klon)
[878]	122	real zalim(klon)
	123	integer lalim(klon)
	124	real norme(klon)
[542]	125	real f(klon), f0(klon)
	126	real zlevinter(klon)
[878]	127	logical therm
[542]	128	logical first
	129	data first /.false./
	130	save first
[987]	131	c$OMP THREADPRIVATE(first)
[542]	132	cRC
	133
	134	character*2 str2
	135	character*10 str10
	136
[1299]	137	character (len=20) :: modname='calcul_sec'
	138	character (len=80) :: abort_message
	139
	140
[940]	141	! LOGICAL vtest(klon),down
[542]	142
	143	EXTERNAL SCOPY
	144
[940]	145	integer ncorrec
[542]	146	save ncorrec
	147	data ncorrec/0/
[987]	148	c$OMP THREADPRIVATE(ncorrec)
[878]	149
[542]	150	c
	151	c-----------------------------------------------------------------------
	152	c initialisation:
	153	c ---------------
	154	c
	155	sorties=.true.
	156	IF(ngrid.NE.klon) THEN
	157	PRINT*
	158	PRINT*,'STOP dans convadj'
	159	PRINT*,'ngrid =',ngrid
	160	PRINT*,'klon =',klon
	161	ENDIF
	162	c
	163	c-----------------------------------------------------------------------
	164	c incrementation eventuelle de tendances precedentes:
	165	c ---------------------------------------------------
	166
[878]	167	c print*,'0 OK convect8'
[542]	168
	169	DO 1010 l=1,nlay
	170	DO 1015 ig=1,ngrid
	171	zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
	172	zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
	173	zu(ig,l)=pu(ig,l)
	174	zv(ig,l)=pv(ig,l)
	175	zo(ig,l)=po(ig,l)
	176	ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
	177	1015 CONTINUE
	178	1010 CONTINUE
	179
[878]	180	c print*,'1 OK convect8'
[542]	181	c --------------------
	182	c
	183	c
	184	c + + + + + + + + + + +
	185	c
	186	c
	187	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	188	c wh,wt,wo ...
	189	c
	190	c + + + + + + + + + + + zh,zu,zv,zo,rho
	191	c
	192	c
	193	c -------------------- zlev(1)
	194	c \\\\\\\\\\\\\\\\\\\\
	195	c
	196	c
	197
	198	c-----------------------------------------------------------------------
	199	c Calcul des altitudes des couches
	200	c-----------------------------------------------------------------------
	201
	202	do l=2,nlay
	203	do ig=1,ngrid
	204	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
	205	enddo
	206	enddo
	207	do ig=1,ngrid
	208	zlev(ig,1)=0.
	209	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
	210	enddo
	211	do l=1,nlay
	212	do ig=1,ngrid
	213	zlay(ig,l)=pphi(ig,l)/RG
	214	enddo
	215	enddo
	216
	217	c print*,'2 OK convect8'
	218	c-----------------------------------------------------------------------
	219	c Calcul des densites
	220	c-----------------------------------------------------------------------
	221
	222	do l=1,nlay
	223	do ig=1,ngrid
	224	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
	225	enddo
	226	enddo
	227
	228	do l=2,nlay
	229	do ig=1,ngrid
	230	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
	231	enddo
	232	enddo
	233
	234	do k=1,nlay
	235	do l=1,nlay+1
	236	do ig=1,ngrid
	237	wa(ig,k,l)=0.
	238	enddo
	239	enddo
	240	enddo
	241
	242	c print*,'3 OK convect8'
	243	c------------------------------------------------------------------
	244	c Calcul de w2, quarre de w a partir de la cape
	245	c a partir de w2, on calcule wa, vitesse de l'ascendance
	246	c
	247	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
	248	c w2 est stoke dans wa
	249	c
	250	c ATTENTION: dans convect8, on n'utilise le calcule des wa
	251	c independants par couches que pour calculer l'entrainement
	252	c a la base et la hauteur max de l'ascendance.
	253	c
	254	c Indicages:
	255	c l'ascendance provenant du niveau k traverse l'interface l avec
	256	c une vitesse wa(k,l).
	257	c
	258	c --------------------
	259	c
	260	c + + + + + + + + + +
	261	c
	262	c wa(k,l) ---- -------------------- l
	263	c /\
	264	c /\|\|\ + + + + + + + + + +
	265	c \|\|
	266	c \|\| --------------------
	267	c \|\|
	268	c \|\| + + + + + + + + + +
	269	c \|\|
	270	c \|\| --------------------
	271	c \|\|__
	272	c \|___ + + + + + + + + + + k
	273	c
	274	c --------------------
	275	c
	276	c
	277	c
	278	c------------------------------------------------------------------
	279
	280	cCR: ponderation entrainement des couches instables
	281	cdef des entr_star tels que entr=f*entr_star
	282	do l=1,klev
	283	do ig=1,ngrid
	284	entr_star(ig,l)=0.
	285	enddo
	286	enddo
	287	c determination de la longueur de la couche d entrainement
	288	do ig=1,ngrid
	289	lentr(ig)=1
	290	enddo
	291
	292	con ne considere que les premieres couches instables
[878]	293	therm=.false.
[542]	294	do k=nlay-2,1,-1
	295	do ig=1,ngrid
	296	if (ztv(ig,k).gt.ztv(ig,k+1).and.
	297	s ztv(ig,k+1).le.ztv(ig,k+2)) then
[878]	298	lentr(ig)=k+1
	299	therm=.true.
[542]	300	endif
	301	enddo
	302	enddo
[878]	303	climitation de la valeur du lentr
	304	c do ig=1,ngrid
	305	c lentr(ig)=min(5,lentr(ig))
	306	c enddo
[542]	307	c determination du lmin: couche d ou provient le thermique
	308	do ig=1,ngrid
	309	lmin(ig)=1
	310	enddo
	311	do ig=1,ngrid
	312	do l=nlay,2,-1
	313	if (ztv(ig,l-1).gt.ztv(ig,l)) then
	314	lmin(ig)=l-1
	315	endif
	316	enddo
	317	enddo
[878]	318	cinitialisations
	319	do ig=1,ngrid
	320	zalim(ig)=0.
	321	norme(ig)=0.
	322	lalim(ig)=1
	323	enddo
	324	do k=1,klev-1
	325	do ig=1,ngrid
	326	zalim(ig)=zalim(ig)+zlev(ig,k)*MAX(0.,(ztv(ig,k)-ztv(ig,k+1))
	327	s /(zlev(ig,k+1)-zlev(ig,k)))
	328	c s *(zlev(ig,k+1)-zlev(ig,k))
	329	norme(ig)=norme(ig)+MAX(0.,(ztv(ig,k)-ztv(ig,k+1))
	330	s /(zlev(ig,k+1)-zlev(ig,k)))
	331	c s *(zlev(ig,k+1)-zlev(ig,k))
	332	enddo
	333	enddo
	334	do ig=1,ngrid
	335	if (norme(ig).gt.1.e-10) then
	336	zalim(ig)=max(10.*zalim(ig)/norme(ig),zlev(ig,2))
	337	c zalim(ig)=min(zalim(ig),zlev(ig,lentr(ig)))
	338	endif
	339	enddo
	340	cdétermination du lalim correspondant
	341	do k=1,klev-1
	342	do ig=1,ngrid
	343	if ((zalim(ig).gt.zlev(ig,k)).and.(zalim(ig).le.zlev(ig,k+1)))
	344	s then
	345	lalim(ig)=k
	346	endif
	347	enddo
	348	enddo
[542]	349	c
	350	c definition de l'entrainement des couches
	351	do l=1,klev-1
	352	do ig=1,ngrid
	353	if (ztv(ig,l).gt.ztv(ig,l+1).and.
[878]	354	s l.ge.lmin(ig).and.l.lt.lentr(ig)) then
	355	entr_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)
	356	c s *(zlev(ig,l+1)-zlev(ig,l))
	357	s *sqrt(zlev(ig,l+1))
	358	cautre def
	359	c entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
	360	c s /zlev(ig,lentr(ig)+2)))**(3./2.)
[542]	361	endif
	362	enddo
	363	enddo
[878]	364	cnouveau test
	365	c if (therm) then
	366	do l=1,klev-1
	367	do ig=1,ngrid
	368	if (ztv(ig,l).gt.ztv(ig,l+1).and.
	369	s l.ge.lmin(ig).and.l.le.lalim(ig)
	370	s .and.zalim(ig).gt.1.e-10) then
	371	c if (l.le.lentr(ig)) then
	372	c entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
	373	c s /zalim(ig)))**(3./2.)
	374	c write(10,*)zlev(ig,l),entr_star(ig,l)
	375	endif
	376	enddo
	377	enddo
	378	c endif
	379	c pas de thermique si couche 1 stable
[542]	380	do ig=1,ngrid
	381	if (lmin(ig).gt.5) then
	382	do l=1,klev
	383	entr_star(ig,l)=0.
	384	enddo
	385	endif
	386	enddo
	387	c calcul de l entrainement total
	388	do ig=1,ngrid
	389	entr_star_tot(ig)=0.
	390	enddo
	391	do ig=1,ngrid
	392	do k=1,klev
	393	entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
	394	enddo
	395	enddo
[878]	396	c Calcul entrainement normalise
	397	do ig=1,ngrid
	398	if (entr_star_tot(ig).gt.1.e-10) then
	399	c do l=1,lentr(ig)
	400	do l=1,klev
	401	cdef possibles pour entr_star: zdthetadz, dthetadz, zdtheta
	402	entr_star(ig,l)=entr_star(ig,l)/entr_star_tot(ig)
	403	enddo
	404	endif
	405	enddo
[542]	406	c
[878]	407	c print*,'fin calcul entr_star'
[542]	408	do k=1,klev
	409	do ig=1,ngrid
	410	ztva(ig,k)=ztv(ig,k)
	411	enddo
	412	enddo
	413	cRC
	414	c print*,'7 OK convect8'
	415	do k=1,klev+1
	416	do ig=1,ngrid
	417	zw2(ig,k)=0.
	418	fmc(ig,k)=0.
	419	cCR
	420	f_star(ig,k)=0.
	421	cRC
	422	larg_cons(ig,k)=0.
	423	larg_detr(ig,k)=0.
	424	wa_moy(ig,k)=0.
	425	enddo
	426	enddo
	427
	428	c print*,'8 OK convect8'
	429	do ig=1,ngrid
	430	linter(ig)=1.
	431	lmaxa(ig)=1
	432	lmix(ig)=1
	433	wmaxa(ig)=0.
	434	enddo
	435
	436	cCR:
	437	do l=1,nlay-2
	438	do ig=1,ngrid
	439	if (ztv(ig,l).gt.ztv(ig,l+1)
	440	s .and.entr_star(ig,l).gt.1.e-10
	441	s .and.zw2(ig,l).lt.1e-10) then
	442	f_star(ig,l+1)=entr_star(ig,l)
	443	ctest:calcul de dteta
	444	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
	445	s *(zlev(ig,l+1)-zlev(ig,l))
	446	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
	447	larg_detr(ig,l)=0.
	448	else if ((zw2(ig,l).ge.1e-10).and.
	449	s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
	450	f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
	451	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
	452	s *ztv(ig,l))/f_star(ig,l+1)
	453	zw2(ig,l+1)=zw2(ig,l)(f_star(ig,l)/f_star(ig,l+1))*2+
	454	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	455	s *(zlev(ig,l+1)-zlev(ig,l))
	456	endif
	457	c determination de zmax continu par interpolation lineaire
	458	if (zw2(ig,l+1).lt.0.) then
	459	ctest
	460	if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
[878]	461	c print*,'pb linter'
[542]	462	endif
	463	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
	464	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	465	zw2(ig,l+1)=0.
	466	lmaxa(ig)=l
	467	else
	468	if (zw2(ig,l+1).lt.0.) then
[878]	469	c print*,'pb1 zw2<0'
[542]	470	endif
	471	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	472	endif
	473	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	474	c lmix est le niveau de la couche ou w (wa_moy) est maximum
	475	lmix(ig)=l+1
	476	wmaxa(ig)=wa_moy(ig,l+1)
	477	endif
	478	enddo
	479	enddo
[878]	480	c print*,'fin calcul zw2'
[542]	481	c
	482	c Calcul de la couche correspondant a la hauteur du thermique
	483	do ig=1,ngrid
	484	lmax(ig)=lentr(ig)
[878]	485	c lmax(ig)=lalim(ig)
[542]	486	enddo
	487	do ig=1,ngrid
	488	do l=nlay,lentr(ig)+1,-1
[878]	489	c do l=nlay,lalim(ig)+1,-1
[542]	490	if (zw2(ig,l).le.1.e-10) then
	491	lmax(ig)=l-1
	492	endif
	493	enddo
	494	enddo
[878]	495	c pas de thermique si couche 1 stable
[542]	496	do ig=1,ngrid
	497	if (lmin(ig).gt.5) then
	498	lmax(ig)=1
	499	lmin(ig)=1
[878]	500	lentr(ig)=1
	501	lalim(ig)=1
[542]	502	endif
	503	enddo
	504	c
	505	c Determination de zw2 max
	506	do ig=1,ngrid
	507	wmax(ig)=0.
	508	enddo
	509
	510	do l=1,nlay
	511	do ig=1,ngrid
	512	if (l.le.lmax(ig)) then
	513	if (zw2(ig,l).lt.0.)then
[878]	514	c print*,'pb2 zw2<0'
[542]	515	endif
	516	zw2(ig,l)=sqrt(zw2(ig,l))
	517	wmax(ig)=max(wmax(ig),zw2(ig,l))
	518	else
	519	zw2(ig,l)=0.
	520	endif
	521	enddo
	522	enddo
	523
	524	c Longueur caracteristique correspondant a la hauteur des thermiques.
	525	do ig=1,ngrid
	526	zmax(ig)=0.
	527	zlevinter(ig)=zlev(ig,1)
	528	enddo
	529	do ig=1,ngrid
	530	c calcul de zlevinter
	531	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
	532	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
	533	s -zlev(ig,lmax(ig)))
	534	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
	535	enddo
[878]	536	do ig=1,ngrid
	537	c write(8,*)zmax(ig),lmax(ig),lentr(ig),lmin(ig)
	538	enddo
[1299]	539	con stope après les calculs de zmax et wmax
[878]	540	RETURN
[542]	541
[878]	542	c print*,'avant fermeture'
[542]	543	c Fermeture,determination de f
[878]	544	cAttention! entrainement normalisé ou pas?
[542]	545	do ig=1,ngrid
	546	entr_star2(ig)=0.
	547	enddo
	548	do ig=1,ngrid
	549	if (entr_star_tot(ig).LT.1.e-10) then
	550	f(ig)=0.
	551	else
	552	do k=lmin(ig),lentr(ig)
[878]	553	c do k=lmin(ig),lalim(ig)
[542]	554	entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
	555	s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
	556	enddo
	557	c Nouvelle fermeture
	558	f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
[878]	559	s *entr_star2(ig))
	560	c s *entr_star_tot(ig)
[542]	561	ctest
	562	c if (first) then
[878]	563	f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
	564	s *wmax(ig))
[542]	565	c endif
	566	endif
[878]	567	f0(ig)=f(ig)
[542]	568	c first=.true.
	569	enddo
[878]	570	c print*,'apres fermeture'
	571	con stoppe après la fermeture
	572	RETURN
[542]	573	c Calcul de l'entrainement
	574	do k=1,klev
	575	do ig=1,ngrid
	576	entr(ig,k)=f(ig)*entr_star(ig,k)
	577	enddo
	578	enddo
[878]	579	con stoppe après le calcul de entr
	580	c RETURN
	581	cCR:test pour entrainer moins que la masse
	582	c do ig=1,ngrid
	583	c do l=1,lentr(ig)
	584	c if ((entr(ig,l)ptimestep).gt.(0.9masse(ig,l))) then
	585	c entr(ig,l+1)=entr(ig,l+1)+entr(ig,l)
	586	c s -0.9*masse(ig,l)/ptimestep
	587	c entr(ig,l)=0.9*masse(ig,l)/ptimestep
	588	c endif
	589	c enddo
	590	c enddo
	591	cCR: fin test
[542]	592	c Calcul des flux
	593	do ig=1,ngrid
	594	do l=1,lmax(ig)-1
	595	fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
	596	enddo
	597	enddo
	598
	599	cRC
	600
	601
	602	c print*,'9 OK convect8'
	603	c print*,'WA1 ',wa_moy
	604
	605	c determination de l'indice du debut de la mixed layer ou w decroit
	606
	607	c calcul de la largeur de chaque ascendance dans le cas conservatif.
	608	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
[878]	609	c d'une couche est égale à la hauteur de la couche alimentante.
[542]	610	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
	611	c de la vitesse d'entrainement horizontal dans la couche alimentante.
	612
	613	do l=2,nlay
	614	do ig=1,ngrid
	615	if (l.le.lmaxa(ig)) then
	616	zw=max(wa_moy(ig,l),1.e-10)
	617	larg_cons(ig,l)=zmax(ig)*r_aspect
	618	s fmc(ig,l)/(rhobarz(ig,l)zw)
	619	endif
	620	enddo
	621	enddo
	622
	623	do l=2,nlay
	624	do ig=1,ngrid
	625	if (l.le.lmaxa(ig)) then
	626	c if (idetr.eq.0) then
	627	c cette option est finalement en dur.
	628	if ((l_mix*zlev(ig,l)).lt.0.)then
[878]	629	c print,'pb l_mixzlev<0'
[542]	630	endif
[878]	631	cCR: test: nouvelle def de lambda
	632	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	633	if (zw2(ig,l).gt.1.e-10) then
	634	larg_detr(ig,l)=sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
	635	else
[542]	636	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
[878]	637	endif
	638	cRC
[542]	639	c else if (idetr.eq.1) then
	640	c larg_detr(ig,l)=larg_cons(ig,l)
	641	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
	642	c else if (idetr.eq.2) then
	643	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	644	c s *sqrt(wa_moy(ig,l))
	645	c else if (idetr.eq.4) then
	646	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	647	c s *wa_moy(ig,l)
	648	c endif
	649	endif
	650	enddo
	651	enddo
	652
	653	c print*,'10 OK convect8'
	654	c print*,'WA2 ',wa_moy
[878]	655	c calcul de la fraction de la maille concernée par l'ascendance en tenant
[542]	656	c compte de l'epluchage du thermique.
	657	c
	658	cCR def de zmix continu (profil parabolique des vitesses)
	659	do ig=1,ngrid
	660	if (lmix(ig).gt.1.) then
	661	c test
	662	if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	663	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
	664	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	665	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
	666	s then
	667	c
	668	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	669	s ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2)
	670	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	671	s ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))
	672	s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	673	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
	674	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	675	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
	676	else
	677	zmix(ig)=zlev(ig,lmix(ig))
[878]	678	c print*,'pb zmix'
[542]	679	endif
	680	else
	681	zmix(ig)=0.
	682	endif
	683	ctest
	684	if ((zmax(ig)-zmix(ig)).lt.0.) then
	685	zmix(ig)=0.99*zmax(ig)
	686	c print*,'pb zmix>zmax'
	687	endif
	688	enddo
	689	c
	690	c calcul du nouveau lmix correspondant
	691	do ig=1,ngrid
	692	do l=1,klev
	693	if (zmix(ig).ge.zlev(ig,l).and.
	694	s zmix(ig).lt.zlev(ig,l+1)) then
	695	lmix(ig)=l
	696	endif
	697	enddo
	698	enddo
	699	c
	700	do l=2,nlay
	701	do ig=1,ngrid
	702	if(larg_cons(ig,l).gt.1.) then
	703	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
	704	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
	705	s /(r_aspect*zmax(ig))
	706	c test
	707	fraca(ig,l)=max(fraca(ig,l),0.)
	708	fraca(ig,l)=min(fraca(ig,l),0.5)
	709	fracd(ig,l)=1.-fraca(ig,l)
	710	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
	711	else
	712	c wa_moy(ig,l)=0.
	713	fraca(ig,l)=0.
	714	fracc(ig,l)=0.
	715	fracd(ig,l)=1.
	716	endif
	717	enddo
	718	enddo
	719	cCR: calcul de fracazmix
	720	do ig=1,ngrid
	721	fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
	722	s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
	723	s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
	724	s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
	725	enddo
	726	c
	727	do l=2,nlay
	728	do ig=1,ngrid
	729	if(larg_cons(ig,l).gt.1.) then
	730	if (l.gt.lmix(ig)) then
	731	ctest
	732	if (zmax(ig)-zmix(ig).lt.1.e-10) then
	733	c print*,'pb xxx'
	734	xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
	735	else
	736	xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
	737	endif
	738	if (idetr.eq.0) then
	739	fraca(ig,l)=fracazmix(ig)
	740	else if (idetr.eq.1) then
	741	fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
	742	else if (idetr.eq.2) then
	743	fraca(ig,l)=fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
	744	else
	745	fraca(ig,l)=fracazmix(ig)xxx(ig,l)*2
	746	endif
	747	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
	748	fraca(ig,l)=max(fraca(ig,l),0.)
	749	fraca(ig,l)=min(fraca(ig,l),0.5)
	750	fracd(ig,l)=1.-fraca(ig,l)
	751	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
	752	endif
	753	endif
	754	enddo
	755	enddo
	756
[878]	757	c print*,'fin calcul fraca'
[542]	758	c print*,'11 OK convect8'
	759	c print*,'Ea3 ',wa_moy
	760	c------------------------------------------------------------------
	761	c Calcul de fracd, wd
	762	c somme wa - wd = 0
	763	c------------------------------------------------------------------
	764
	765
	766	do ig=1,ngrid
	767	fm(ig,1)=0.
	768	fm(ig,nlay+1)=0.
	769	enddo
	770
	771	do l=2,nlay
	772	do ig=1,ngrid
	773	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
	774	cCR:test
	775	if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
	776	s .and.l.gt.lmix(ig)) then
	777	fm(ig,l)=fm(ig,l-1)
	778	c write(1,*)'ajustement fm, l',l
	779	endif
	780	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
	781	cRC
	782	enddo
	783	do ig=1,ngrid
	784	if(fracd(ig,l).lt.0.1) then
[1299]	785	abort_message = 'fracd trop petit'
	786	CALL abort_gcm (modname,abort_message,1)
	787
[542]	788	else
	789	c vitesse descendante "diagnostique"
	790	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
	791	endif
	792	enddo
	793	enddo
	794
	795	do l=1,nlay
	796	do ig=1,ngrid
	797	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
	798	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
	799	enddo
	800	enddo
	801
[878]	802	c print*,'12 OK convect8'
[542]	803	c print*,'WA4 ',wa_moy
	804	cc------------------------------------------------------------------
	805	c calcul du transport vertical
	806	c------------------------------------------------------------------
	807
	808	go to 4444
	809	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
	810	do l=2,nlay-1
	811	do ig=1,ngrid
	812	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
	813	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
	814	c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
	815	c s ,fm(ig,l+1)*ptimestep
	816	c s ,' M=',masse(ig,l),masse(ig,l+1)
	817	endif
	818	enddo
	819	enddo
	820
	821	do l=1,nlay
	822	do ig=1,ngrid
	823	if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
	824	c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
	825	c s ,entr(ig,l)*ptimestep
	826	c s ,' M=',masse(ig,l)
	827	endif
	828	enddo
	829	enddo
	830
	831	do l=1,nlay
	832	do ig=1,ngrid
	833	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
	834	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
	835	c s ,' FM=',fm(ig,l)
	836	endif
	837	if(.not.masse(ig,l).ge.1.e-10
	838	s .or..not.masse(ig,l).le.1.e4) then
	839	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
	840	c s ,' M=',masse(ig,l)
	841	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
	842	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
	843	c print*,'zlev(ig,l+1),zlev(ig,l)'
	844	c s ,zlev(ig,l+1),zlev(ig,l)
	845	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
	846	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
	847	endif
	848	if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
	849	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
	850	c s ,' E=',entr(ig,l)
	851	endif
	852	enddo
	853	enddo
	854
	855	4444 continue
	856
	857	cCR:redefinition du entr
	858	do l=1,nlay
	859	do ig=1,ngrid
	860	detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
	861	if (detr(ig,l).lt.0.) then
[878]	862	c entr(ig,l)=entr(ig,l)-detr(ig,l)
	863	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
[542]	864	detr(ig,l)=0.
	865	c print*,'WARNING !!! detrainement negatif ',ig,l
	866	endif
	867	enddo
	868	enddo
	869	cRC
	870	if (w2di.eq.1) then
[1299]	871	fm0=fm0+ptimestep*(fm-fm0)/REAL(tho)
	872	entr0=entr0+ptimestep*(entr-entr0)/REAL(tho)
[542]	873	else
	874	fm0=fm
	875	entr0=entr
	876	endif
	877
	878	if (1.eq.1) then
	879	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	880	. ,zh,zdhadj,zha)
	881	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	882	. ,zo,pdoadj,zoa)
	883	else
	884	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
	885	. ,zh,zdhadj,zha)
	886	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
	887	. ,zo,pdoadj,zoa)
	888	endif
	889
	890	if (1.eq.0) then
	891	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
	892	. ,fraca,zmax
	893	. ,zu,zv,pduadj,pdvadj,zua,zva)
	894	else
	895	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	896	. ,zu,pduadj,zua)
	897	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	898	. ,zv,pdvadj,zva)
	899	endif
	900
	901	do l=1,nlay
	902	do ig=1,ngrid
	903	zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
	904	zf2=zf/(1.-zf)
	905	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
	906	wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
	907	enddo
	908	enddo
	909
	910
	911
	912	c print*,'13 OK convect8'
	913	c print*,'WA5 ',wa_moy
	914	do l=1,nlay
	915	do ig=1,ngrid
	916	pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
	917	enddo
	918	enddo
	919
	920
	921	c do l=1,nlay
	922	c do ig=1,ngrid
	923	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
	924	c print*,'WARN!!! ig=',ig,' l=',l
	925	c s ,' pdtadj=',pdtadj(ig,l)
	926	c endif
	927	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
	928	c print*,'WARN!!! ig=',ig,' l=',l
	929	c s ,' pdoadj=',pdoadj(ig,l)
	930	c endif
	931	c enddo
	932	c enddo
	933
[878]	934	c print*,'14 OK convect8'
[542]	935	c------------------------------------------------------------------
	936	c Calculs pour les sorties
	937	c------------------------------------------------------------------
	938
	939	if(sorties) then
	940	do l=1,nlay
	941	do ig=1,ngrid
	942	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
	943	zld(ig,l)=fracd(ig,l)*zmax(ig)
	944	if(1.-fracd(ig,l).gt.1.e-10)
	945	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
	946	enddo
	947	enddo
	948
	949	cdeja fait
	950	c do l=1,nlay
	951	c do ig=1,ngrid
	952	c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
	953	c if (detr(ig,l).lt.0.) then
	954	c entr(ig,l)=entr(ig,l)-detr(ig,l)
	955	c detr(ig,l)=0.
	956	c print*,'WARNING !!! detrainement negatif ',ig,l
	957	c endif
	958	c enddo
	959	c enddo
	960
	961	c print*,'15 OK convect8'
	962
	963	isplit=isplit+1
	964
	965
	966	c #define und
	967	goto 123
	968	#ifdef und
	969	CALL writeg1d(1,nlay,wd,'wd ','wd ')
	970	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
	971	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
	972	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
	973	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
	974	CALL writeg1d(1,nlay,zla,'la ','la ')
	975	CALL writeg1d(1,nlay,zld,'ld ','ld ')
	976	CALL writeg1d(1,nlay,pt,'pt ','pt ')
	977	CALL writeg1d(1,nlay,zh,'zh ','zh ')
	978	CALL writeg1d(1,nlay,zha,'zha ','zha ')
	979	CALL writeg1d(1,nlay,zu,'zu ','zu ')
	980	CALL writeg1d(1,nlay,zv,'zv ','zv ')
	981	CALL writeg1d(1,nlay,zo,'zo ','zo ')
	982	CALL writeg1d(1,nlay,wh,'wh ','wh ')
	983	CALL writeg1d(1,nlay,wu,'wu ','wu ')
	984	CALL writeg1d(1,nlay,wv,'wv ','wv ')
	985	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
	986	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
	987	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
	988	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
	989	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
	990	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
	991	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
	992	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
	993
	994	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
	995	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
	996	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
	997
	998	c recalcul des flux en diagnostique...
	999	c print*,'PAS DE TEMPS ',ptimestep
	1000	call dt2F(pplev,pplay,pt,pdtadj,wh)
	1001	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
	1002	#endif
	1003	123 continue
[878]	1004	#define troisD
[542]	1005	#ifdef troisD
	1006	c if (sorties) then
[938]	1007	c print*,'Debut des wrgradsfi'
[542]	1008
	1009	c print*,'16 OK convect8'
	1010	call wrgradsfi(1,nlay,wd,'wd ','wd ')
	1011	call wrgradsfi(1,nlay,zwa,'wa ','wa ')
	1012	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
	1013	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
	1014	call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
	1015	call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
	1016	c print*,'WA6 ',wa_moy
	1017	call wrgradsfi(1,nlay,zla,'la ','la ')
	1018	call wrgradsfi(1,nlay,zld,'ld ','ld ')
	1019	call wrgradsfi(1,nlay,pt,'pt ','pt ')
	1020	call wrgradsfi(1,nlay,zh,'zh ','zh ')
	1021	call wrgradsfi(1,nlay,zha,'zha ','zha ')
	1022	call wrgradsfi(1,nlay,zua,'zua ','zua ')
	1023	call wrgradsfi(1,nlay,zva,'zva ','zva ')
	1024	call wrgradsfi(1,nlay,zu,'zu ','zu ')
	1025	call wrgradsfi(1,nlay,zv,'zv ','zv ')
	1026	call wrgradsfi(1,nlay,zo,'zo ','zo ')
	1027	call wrgradsfi(1,nlay,wh,'wh ','wh ')
	1028	call wrgradsfi(1,nlay,wu,'wu ','wu ')
	1029	call wrgradsfi(1,nlay,wv,'wv ','wv ')
	1030	call wrgradsfi(1,nlay,wo,'wo ','wo ')
	1031	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
	1032	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
	1033	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
	1034	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
	1035	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
	1036	call wrgradsfi(1,nlay,entr,'entr ','entr ')
	1037	call wrgradsfi(1,nlay,detr,'detr ','detr ')
	1038	call wrgradsfi(1,nlay,fm,'fm ','fm ')
	1039	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
	1040	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
	1041	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
	1042	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
	1043
	1044	call wrgradsfi(1,nlay,zo,'zo ','zo ')
	1045	call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
	1046	call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
	1047
	1048	cCR:nouveaux diagnostiques
	1049	call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')
	1050	call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')
	1051	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
	1052	call wrgradsfi(1,1,zmix,'zmix ','zmix ')
	1053	zsortie1d(:)=lmax(:)
	1054	call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')
	1055	call wrgradsfi(1,1,wmax,'wmax ','wmax ')
	1056	zsortie1d(:)=lmix(:)
	1057	call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')
	1058	zsortie1d(:)=lentr(:)
	1059	call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')
	1060
	1061	c print*,'17 OK convect8'
	1062
	1063	do k=1,klev/10
	1064	write(str2,'(i2.2)') k
	1065	str10='wa'//str2
	1066	do l=1,nlay
	1067	do ig=1,ngrid
	1068	zsortie(ig,l)=wa(ig,k,l)
	1069	enddo
	1070	enddo
	1071	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
	1072	do l=1,nlay
	1073	do ig=1,ngrid
	1074	zsortie(ig,l)=larg_part(ig,k,l)
	1075	enddo
	1076	enddo
	1077	str10='la'//str2
	1078	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
	1079	enddo
	1080
	1081
	1082	c print*,'18 OK convect8'
	1083	c endif
[938]	1084	c print*,'Fin des wrgradsfi'
[542]	1085	#endif
	1086
	1087	endif
	1088
	1089	c if(wa_moy(1,4).gt.1.e-10) stop
	1090
[878]	1091	c print*,'19 OK convect8'
[542]	1092	return
	1093	end
	1094
[878]	1095	SUBROUTINE fermeture_seche(ngrid,nlay
	1096	s ,pplay,pplev,pphi,zlev,rhobarz,f0,zpspsk
	1097	s ,alim_star,zh,zo,lentr,lmin,nu_min,nu_max,r_aspect
	1098	s ,zmax,wmax)
[542]	1099
[940]	1100	USE dimphy
[878]	1101	IMPLICIT NONE
[542]	1102
[878]	1103	#include "dimensions.h"
[940]	1104	cccc#include "dimphy.h"
[878]	1105	#include "YOMCST.h"
[542]	1106
[878]	1107	INTEGER ngrid,nlay
	1108	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	1109	real pphi(ngrid,nlay)
	1110	real zlev(klon,klev+1)
	1111	real alim_star(klon,klev)
	1112	real f0(klon)
	1113	integer lentr(klon)
	1114	integer lmin(klon)
	1115	real zmax(klon)
	1116	real wmax(klon)
	1117	real nu_min
	1118	real nu_max
	1119	real r_aspect
	1120	real rhobarz(klon,klev+1)
	1121	REAL zh(klon,klev)
	1122	real zo(klon,klev)
	1123	real zpspsk(klon,klev)
[542]	1124
[878]	1125	integer ig,l
[542]	1126
[878]	1127	real f_star(klon,klev+1)
	1128	real detr_star(klon,klev)
	1129	real entr_star(klon,klev)
	1130	real zw2(klon,klev+1)
	1131	real linter(klon)
	1132	integer lmix(klon)
	1133	integer lmax(klon)
	1134	real zlevinter(klon)
	1135	real wa_moy(klon,klev+1)
	1136	real wmaxa(klon)
	1137	REAL ztv(klon,klev)
	1138	REAL ztva(klon,klev)
	1139	real nu(klon,klev)
[940]	1140	! real zmax0_sec(klon)
	1141	! save zmax0_sec
	1142	REAL, SAVE, ALLOCATABLE :: zmax0_sec(:)
	1143	c$OMP THREADPRIVATE(zmax0_sec)
	1144	logical, save :: first = .true.
[987]	1145	c$OMP THREADPRIVATE(first)
[542]	1146
[940]	1147	if (first) then
	1148	allocate(zmax0_sec(klon))
	1149	first=.false.
	1150	endif
	1151
[878]	1152	do l=1,nlay
[542]	1153	do ig=1,ngrid
[878]	1154	ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
	1155	ztv(ig,l)=ztv(ig,l)(1.+RETVzo(ig,l))
[542]	1156	enddo
	1157	enddo
[878]	1158	do l=1,nlay-2
[542]	1159	do ig=1,ngrid
[878]	1160	if (ztv(ig,l).gt.ztv(ig,l+1)
	1161	s .and.alim_star(ig,l).gt.1.e-10
	1162	s .and.zw2(ig,l).lt.1e-10) then
	1163	f_star(ig,l+1)=alim_star(ig,l)
	1164	ctest:calcul de dteta
	1165	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
	1166	s *(zlev(ig,l+1)-zlev(ig,l))
	1167	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
	1168	else if ((zw2(ig,l).ge.1e-10).and.
	1169	s (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
	1170	cestimation du detrainement a partir de la geometrie du pas precedent
	1171	ctests sur la definition du detr
	1172	nu(ig,l)=(nu_min+nu_max)/2.
	1173	s *(1.-(nu_max-nu_min)/(nu_max+nu_min)
	1174	s *tanh((((ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l))/0.0005)))
	1175
	1176	detr_star(ig,l)=rhobarz(ig,l)
	1177	s *sqrt(zw2(ig,l))
	1178	s /(r_aspectzmax0_sec(ig))
	1179	c s /(r_aspectzmax0(ig))
	1180	s (sqrt(nu(ig,l)*zlev(ig,l+1)
	1181	s /sqrt(zw2(ig,l)))
	1182	s -sqrt(nu(ig,l)*zlev(ig,l)
	1183	s /sqrt(zw2(ig,l))))
	1184	detr_star(ig,l)=detr_star(ig,l)/f0(ig)
	1185	if ((detr_star(ig,l)).gt.f_star(ig,l)) then
	1186	detr_star(ig,l)=f_star(ig,l)
	1187	endif
	1188	entr_star(ig,l)=0.9*detr_star(ig,l)
	1189	if ((l.lt.lentr(ig))) then
	1190	entr_star(ig,l)=0.
	1191	c detr_star(ig,l)=0.
	1192	endif
[938]	1193	c print*,'ok detr_star'
[878]	1194	cprise en compte du detrainement dans le calcul du flux
	1195	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
	1196	s -detr_star(ig,l)
	1197	ctest sur le signe de f_star
	1198	if ((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10) then
	1199	cAM on melange Tl et qt du thermique
	1200	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+(entr_star(ig,l)
	1201	s +alim_star(ig,l))
	1202	s *ztv(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
	1203	zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)
	1204	s /(f_star(ig,l+1)+detr_star(ig,l)))**2+
	1205	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	1206	s *(zlev(ig,l+1)-zlev(ig,l))
[542]	1207	endif
[878]	1208	endif
[542]	1209	c
[878]	1210	if (zw2(ig,l+1).lt.0.) then
	1211	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
	1212	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	1213	zw2(ig,l+1)=0.
[938]	1214	c print*,'linter=',linter(ig)
[542]	1215	else
[878]	1216	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
[542]	1217	endif
[878]	1218	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	1219	c lmix est le niveau de la couche ou w (wa_moy) est maximum
	1220	lmix(ig)=l+1
	1221	wmaxa(ig)=wa_moy(ig,l+1)
	1222	endif
[542]	1223	enddo
	1224	enddo
[938]	1225	c print*,'fin calcul zw2'
[878]	1226	c
	1227	c Calcul de la couche correspondant a la hauteur du thermique
[542]	1228	do ig=1,ngrid
[878]	1229	lmax(ig)=lentr(ig)
[542]	1230	enddo
	1231	do ig=1,ngrid
[878]	1232	do l=nlay,lentr(ig)+1,-1
	1233	if (zw2(ig,l).le.1.e-10) then
	1234	lmax(ig)=l-1
[542]	1235	endif
	1236	enddo
	1237	enddo
[878]	1238	c pas de thermique si couche 1 stable
[542]	1239	do ig=1,ngrid
[878]	1240	if (lmin(ig).gt.1) then
	1241	lmax(ig)=1
	1242	lmin(ig)=1
	1243	lentr(ig)=1
	1244	endif
	1245	enddo
	1246	c
	1247	c Determination de zw2 max
[542]	1248	do ig=1,ngrid
[878]	1249	wmax(ig)=0.
[542]	1250	enddo
	1251
[878]	1252	do l=1,nlay
[542]	1253	do ig=1,ngrid
[878]	1254	if (l.le.lmax(ig)) then
	1255	if (zw2(ig,l).lt.0.)then
[938]	1256	c print*,'pb2 zw2<0'
[878]	1257	endif
	1258	zw2(ig,l)=sqrt(zw2(ig,l))
	1259	wmax(ig)=max(wmax(ig),zw2(ig,l))
[542]	1260	else
[878]	1261	zw2(ig,l)=0.
[542]	1262	endif
[878]	1263	enddo
[542]	1264	enddo
	1265
[878]	1266	c Longueur caracteristique correspondant a la hauteur des thermiques.
	1267	do ig=1,ngrid
	1268	zmax(ig)=0.
	1269	zlevinter(ig)=zlev(ig,1)
[542]	1270	enddo
[878]	1271	do ig=1,ngrid
	1272	c calcul de zlevinter
	1273	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
	1274	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
	1275	s -zlev(ig,lmax(ig)))
	1276	cpour le cas ou on prend tjs lmin=1
	1277	c zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
	1278	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
	1279	zmax0_sec(ig)=zmax(ig)
[542]	1280	enddo
[878]	1281
	1282	RETURN
	1283	END

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

Download in other formats:

Original Format