Context Navigation

calltherm.F90 @ 1320

Last change on this file since 1320 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 11.7 KB

Rev	Line
[878]	1	!
	2	! $Header$
	3	!
	4	subroutine calltherm(dtime &
	5	& ,pplay,paprs,pphi,weak_inversion &
	6	& ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
	7	& ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
[973]	8	& ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth, &
[927]	9	& ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
[1026]	10	& zmax0,f0,zw2,fraca)
[878]	11
[940]	12	USE dimphy
[878]	13	implicit none
	14	#include "dimensions.h"
[940]	15	!#include "dimphy.h"
[878]	16	#include "thermcell.h"
[938]	17	#include "iniprint.h"
[878]	18
[973]	19	!IM 140508
	20	INTEGER itap
[878]	21	REAL dtime
	22	LOGICAL debut
[973]	23	LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
	24	REAL fact(klon)
	25	INTEGER nbptspb
	26
[878]	27	REAL u_seri(klon,klev),v_seri(klon,klev)
	28	REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
	29	REAL weak_inversion(klon)
	30	REAL paprs(klon,klev+1)
	31	REAL pplay(klon,klev)
	32	REAL pphi(klon,klev)
	33	real zlev(klon,klev+1)
[879]	34	!test: on sort lentr et a* pour alimenter KE
	35	REAL wght_th(klon,klev)
	36	INTEGER lalim_conv(klon)
[1026]	37	REAL zw2(klon,klev+1),fraca(klon,klev+1)
[878]	38
	39	!FH Update Thermiques
	40	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
	41	REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
[973]	42	real fm_therm(klon,klev+1)
	43	real entr_therm(klon,klev),detr_therm(klon,klev)
[878]	44
	45	!********************************************************
	46	! declarations
[973]	47	real fmc_therm(klon,klev+1),zqasc(klon,klev)
[878]	48	real zqla(klon,klev)
[1026]	49	real zqta(klon,klev)
[878]	50	real wmax_sec(klon)
	51	real zmax_sec(klon)
	52	real f_sec(klon)
[973]	53	real detrc_therm(klon,klev)
	54	! FH WARNING : il semble que ces save ne servent a rien
	55	! save fmc_therm, detrc_therm
[878]	56	real clwcon0(klon,klev)
	57	real zqsat(klon,klev)
	58	real zw_sec(klon,klev+1)
	59	integer lmix_sec(klon)
	60	integer lmax(klon)
	61	real ratqscth(klon,klev)
	62	real ratqsdiff(klon,klev)
	63	real zqsatth(klon,klev)
[879]	64	!nouvelles variables pour la convection
	65	real Ale_bl(klon)
	66	real Alp_bl(klon)
	67	real Ale(klon)
	68	real Alp(klon)
	69	!RC
[927]	70	!on garde le zmax du pas de temps precedent
	71	real zmax0(klon), f0(klon)
[878]	72	!********************************************************
	73
	74
	75	! variables locales
	76	REAL d_t_the(klon,klev), d_q_the(klon,klev)
	77	REAL d_u_the(klon,klev),d_v_the(klon,klev)
	78	!
[973]	79	real zfm_therm(klon,klev+1),zdt
	80	real zentr_therm(klon,klev),zdetr_therm(klon,klev)
	81	! FH A VERIFIER : SAVE INUTILES
[940]	82	! save zentr_therm,zfm_therm
[973]	83
[878]	84	integer i,k
[973]	85	logical, save :: first=.true.
[987]	86	!$OMP THREADPRIVATE(first)
[878]	87	!********************************************************
[973]	88	if (first) then
	89	itap=0
	90	first=.false.
	91	endif
[878]	92
[973]	93	! Incrementer le compteur de la physique
	94	itap = itap + 1
	95
[878]	96	! Modele du thermique
	97	! ===================
	98	! print*,'thermiques: WARNING on passe t au lieu de t_seri'
	99
	100
[973]	101	! On prend comme valeur initiale des thermiques la valeur du pas
	102	! de temps precedent
	103	zfm_therm(:,:)=fm_therm(:,:)
	104	zdetr_therm(:,:)=detr_therm(:,:)
	105	zentr_therm(:,:)=entr_therm(:,:)
	106
	107	! On reinitialise les flux de masse a zero pour le cumul en
	108	! cas de splitting
[878]	109	fm_therm(:,:)=0.
	110	entr_therm(:,:)=0.
[973]	111	detr_therm(:,:)=0.
	112
[879]	113	Ale_bl(:)=0.
	114	Alp_bl(:)=0.
[938]	115	if (prt_level.ge.10) then
	116	print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
	117	endif
[878]	118
	119	! tests sur les valeurs negatives de l'eau
[973]	120	logexpr0=prt_level.ge.10
	121	nbptspb=0
[878]	122	do k=1,klev
	123	do i=1,klon
[1146]	124	! Attention teste abderr 19-03-09
	125	! logexpr2(i,k)=.not.q_seri(i,k).ge.0.
	126	logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
[973]	127	if (logexpr2(i,k)) then
[938]	128	q_seri(i,k)=1.e-15
[973]	129	nbptspb=nbptspb+1
[878]	130	endif
[973]	131	! if (logexpr0) &
	132	! & print*,'WARN eau<0 avant therm i=',i,' k=',k &
	133	! & ,' dq,q',d_q_the(i,k),q_seri(i,k)
[878]	134	enddo
	135	enddo
[973]	136	if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
[878]	137
	138	zdt=dtime/float(nsplit_thermals)
	139	do isplit=1,nsplit_thermals
	140
	141	if (iflag_thermals.eq.1) then
	142	CALL thermcell_2002(klon,klev,zdt &
	143	& ,pplay,paprs,pphi &
	144	& ,u_seri,v_seri,t_seri,q_seri &
	145	& ,d_u_the,d_v_the,d_t_the,d_q_the &
	146	& ,zfm_therm,zentr_therm &
	147	& ,r_aspect_thermals,30.,w2di_thermals &
[973]	148	& ,tau_thermals,3)
[878]	149	else if (iflag_thermals.eq.2) then
	150	CALL thermcell_sec(klon,klev,zdt &
	151	& ,pplay,paprs,pphi,zlev &
	152	& ,u_seri,v_seri,t_seri,q_seri &
	153	& ,d_u_the,d_v_the,d_t_the,d_q_the &
	154	& ,zfm_therm,zentr_therm &
	155	& ,r_aspect_thermals,30.,w2di_thermals &
[973]	156	& ,tau_thermals,3)
[878]	157	else if (iflag_thermals.eq.3) then
	158	CALL thermcell(klon,klev,zdt &
	159	& ,pplay,paprs,pphi &
	160	& ,u_seri,v_seri,t_seri,q_seri &
	161	& ,d_u_the,d_v_the,d_t_the,d_q_the &
	162	& ,zfm_therm,zentr_therm &
	163	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
[973]	164	& ,tau_thermals,3)
[878]	165	else if (iflag_thermals.eq.10) then
	166	CALL thermcell_eau(klon,klev,zdt &
	167	& ,pplay,paprs,pphi &
	168	& ,u_seri,v_seri,t_seri,q_seri &
	169	& ,d_u_the,d_v_the,d_t_the,d_q_the &
	170	& ,zfm_therm,zentr_therm &
	171	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
[973]	172	& ,tau_thermals,3)
[878]	173	else if (iflag_thermals.eq.11) then
[1146]	174	stop 'cas non prevu dans calltherm'
[878]	175	! CALL thermcell_pluie(klon,klev,zdt &
	176	! & ,pplay,paprs,pphi,zlev &
	177	! & ,u_seri,v_seri,t_seri,q_seri &
	178	! & ,d_u_the,d_v_the,d_t_the,d_q_the &
	179	! & ,zfm_therm,zentr_therm,zqla &
	180	! & ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
[973]	181	! & ,tau_thermals,3)
[878]	182	else if (iflag_thermals.eq.12) then
	183	CALL calcul_sec(klon,klev,zdt &
	184	& ,pplay,paprs,pphi,zlev &
	185	& ,u_seri,v_seri,t_seri,q_seri &
	186	& ,zmax_sec,wmax_sec,zw_sec,lmix_sec &
	187	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
[973]	188	& ,tau_thermals)
[878]	189	else if (iflag_thermals.ge.13) then
[973]	190	CALL thermcell_main(itap,klon,klev,zdt &
[878]	191	& ,pplay,paprs,pphi,debut &
	192	& ,u_seri,v_seri,t_seri,q_seri &
	193	& ,d_u_the,d_v_the,d_t_the,d_q_the &
[1026]	194	& ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax &
[878]	195	& ,ratqscth,ratqsdiff,zqsatth &
[973]	196	& ,r_aspect_thermals,l_mix_thermals &
	197	& ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
[1026]	198	& ,zmax0,f0,zw2,fraca)
[878]	199	endif
	200
	201
[973]	202	fact(:)=0.
[878]	203	DO i=1,klon
[973]	204	logexpr1(i)=iflag_thermals.lt.14.or.weak_inversion(i).gt.0.5
	205	IF(logexpr1(i)) fact(i)=1./float(nsplit_thermals)
	206	ENDDO
[878]	207
[973]	208	DO k=1,klev
[878]	209	! transformation de la derivee en tendance
[973]	210	d_t_the(:,k)=d_t_the(:,k)dtimefact(:)
	211	d_u_the(:,k)=d_u_the(:,k)dtimefact(:)
	212	d_v_the(:,k)=d_v_the(:,k)dtimefact(:)
	213	d_q_the(:,k)=d_q_the(:,k)dtimefact(:)
	214	fm_therm(:,k)=fm_therm(:,k) &
	215	& +zfm_therm(:,k)*fact(:)
	216	entr_therm(:,k)=entr_therm(:,k) &
	217	& +zentr_therm(:,k)*fact(:)
[1026]	218	detr_therm(:,k)=detr_therm(:,k) &
	219	& +zdetr_therm(:,k)*fact(:)
[973]	220	ENDDO
	221	fm_therm(:,klev+1)=0.
[878]	222
	223
	224
	225	! accumulation de la tendance
[973]	226	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
	227	d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
	228	d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
	229	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
[878]	230
	231	! incrementation des variables meteo
[973]	232	t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
	233	u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
	234	v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
	235	qmemoire(:,:)=q_seri(:,:)
	236	q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
[878]	237
[879]	238	DO i=1,klon
[973]	239	if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
[879]	240	fm_therm(i,klev+1)=0.
	241	Ale_bl(i)=Ale_bl(i)+Ale(i)/float(nsplit_thermals)
	242	! write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
	243	Alp_bl(i)=Alp_bl(i)+Alp(i)/float(nsplit_thermals)
	244	! write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
	245	ENDDO
	246
[973]	247	!IM 060508 marche pas comme cela !!! enddo ! isplit
	248
[878]	249	! tests sur les valeurs negatives de l'eau
[973]	250	nbptspb=0
[878]	251	DO k = 1, klev
	252	DO i = 1, klon
[973]	253	logexpr2(i,k)=.not.q_seri(i,k).ge.0.
	254	if (logexpr2(i,k)) then
	255	q_seri(i,k)=1.e-15
	256	nbptspb=nbptspb+1
	257	! if (prt_level.ge.10) then
	258	! print*,'WARN eau<0 apres therm i=',i,' k=',k &
	259	! & ,' dq,q',d_q_the(i,k),q_seri(i,k), &
	260	! & 'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
[938]	261	endif
[878]	262	! stop
	263	ENDDO
	264	ENDDO
[973]	265	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
[878]	266	! tests sur les valeurs de la temperature
[973]	267	nbptspb=0
[878]	268	DO k = 1, klev
	269	DO i = 1, klon
[973]	270	logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
	271	if (logexpr2(i,k)) nbptspb=nbptspb+1
	272	! if ((t_seri(i,k).lt.50.) .or. &
	273	! & (t_seri(i,k).gt.370.)) then
	274	! print*,'WARN temp apres therm i=',i,' k=',k &
	275	! & ,' t_seri',t_seri(i,k)
[878]	276	! CALL abort
[973]	277	! endif
[878]	278	ENDDO
	279	ENDDO
[973]	280	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
[878]	281	enddo ! isplit
	282
	283	!
	284	!***************************************************************
	285	! calcul du flux ascencant conservatif
	286	! print*,'<<<<calcul flux ascendant conservatif'
	287
	288	fmc_therm=0.
	289	do k=1,klev
	290	do i=1,klon
	291	if (entr_therm(i,k).gt.0.) then
	292	fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
	293	else
	294	fmc_therm(i,k+1)=fmc_therm(i,k)
	295	endif
	296	detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1)) &
	297	& -(fmc_therm(i,k)-fm_therm(i,k))
	298	enddo
	299	enddo
	300
	301
	302	!****************************************************************
	303	! calcul de l'humidite dans l'ascendance
	304	! print*,'<<<<calcul de lhumidite dans thermique'
	305	!CR:on ne le calcule que pour le cas sec
	306	if (iflag_thermals.le.11) then
	307	do i=1,klon
	308	zqasc(i,1)=q_seri(i,1)
	309	do k=2,klev
	310	if (fmc_therm(i,k+1).gt.1.e-6) then
	311	zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1) &
	312	& +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
	313	!CR:test on asseche le thermique
	314	! zqasc(i,k)=zqasc(i,k)/2.
	315	! else
	316	! zqasc(i,k)=q_seri(i,k)
	317	endif
	318	enddo
	319	enddo
	320
	321
	322	! calcul de l'eau condensee dans l'ascendance
	323	! print*,'<<<<calcul de leau condensee dans thermique'
	324	do i=1,klon
	325	do k=1,klev
	326	clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
	327	if (clwcon0(i,k).lt.0. .or. &
	328	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
	329	clwcon0(i,k)=0.
	330	endif
	331	enddo
	332	enddo
	333	else
	334	do i=1,klon
	335	do k=1,klev
	336	clwcon0(i,k)=zqla(i,k)
	337	if (clwcon0(i,k).lt.0. .or. &
	338	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
	339	clwcon0(i,k)=0.
	340	endif
	341	enddo
	342	enddo
	343	endif
	344	!*******************************************************************
	345
	346
	347	return
	348
	349	end

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

Download in other formats:

Original Format