Context Navigation

calltherm.F90 @ 5452

Last change on this file since 5452 was 1428, checked in by Laurent Fairhead, 14 years ago

Problems with some initialisations in the new physics package. These are
temporary fixes. Further investigations needed. Trac ticket issued. JYG

Problèmes d'initialisations dans la nouvelle physique. Ces corrections
sont temporaires, il faut trouver la vraie raison des plantages. Un ticket
trac est créé. JYG

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

File size: 13.2 KB

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

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

Download in other formats:

Original Format