Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

yamada4.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 22.7 KB

Rev	Line
[1]	1	!
	2	! $Header$
	3	!
	4	SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp &
	5	& ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar &
	6	& ,iflag_pbl)
	7	use dimphy
	8	IMPLICIT NONE
	9	#include "iniprint.h"
	10	!c.......................................................................
	11	!cym#include "dimensions.h"
	12	!cym#include "dimphy.h"
	13	!c.......................................................................
	14	!c
	15	!c dt : pas de temps
	16	!c g : g
	17	!c zlev : altitude a chaque niveau (interface inferieure de la couche
	18	!c de meme indice)
	19	!c zlay : altitude au centre de chaque couche
	20	!c u,v : vitesse au centre de chaque couche
	21	!c (en entree : la valeur au debut du pas de temps)
	22	!c teta : temperature potentielle au centre de chaque couche
	23	!c (en entree : la valeur au debut du pas de temps)
	24	!c cd : cdrag
	25	!c (en entree : la valeur au debut du pas de temps)
	26	!c q2 : $q^2$ au bas de chaque couche
	27	!c (en entree : la valeur au debut du pas de temps)
	28	!c (en sortie : la valeur a la fin du pas de temps)
	29	!c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
	30	!c couche)
	31	!c (en sortie : la valeur a la fin du pas de temps)
	32	!c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
	33	!c (en sortie : la valeur a la fin du pas de temps)
	34	!c
	35	!c iflag_pbl doit valoir entre 6 et 9
	36	!c l=6, on prend systematiquement une longueur d'equilibre
	37	!c iflag_pbl=6 : MY 2.0
	38	!c iflag_pbl=7 : MY 2.0.Fournier
	39	!c iflag_pbl=8/9 : MY 2.5
	40	!c iflag_pbl=8 with special obsolete treatments for convergence
	41	!c with Cmpi5 NPv3.1 simulations
	42	!c iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact
	43	!c iflag_pbl=12 = 11 with vertical diffusion off q2
	44	!c
	45	!c 2013/04/01 (FH hourdin@lmd.jussieu.fr)
	46	!c Correction for very stable PBLs (iflag_pbl=10 and 11)
	47	!c iflag_pbl=8 converges numerically with NPv3.1
	48	!c iflag_pbl=11 -> the model starts with NP from start files created by ce0l
	49	!c -> the model can run with longer time-steps.
	50	!c.......................................................................
	51
	52	REAL dt,g,rconst
	53	real plev(klon,klev+1),temp(klon,klev)
	54	real ustar(klon)
	55	real kmin,qmin,pblhmin(klon),coriol(klon)
	56	REAL zlev(klon,klev+1)
	57	REAL zlay(klon,klev)
	58	REAL u(klon,klev)
	59	REAL v(klon,klev)
	60	REAL teta(klon,klev)
	61	REAL cd(klon)
	62	REAL q2(klon,klev+1),qpre
	63	REAL unsdz(klon,klev)
	64	REAL unsdzdec(klon,klev+1)
	65
	66	REAL km(klon,klev+1)
	67	REAL kmpre(klon,klev+1),tmp2
	68	REAL mpre(klon,klev+1)
	69	REAL kn(klon,klev+1)
	70	REAL kq(klon,klev+1)
	71	real ff(klon,klev+1),delta(klon,klev+1)
	72	real aa(klon,klev+1),aa0,aa1
	73	integer iflag_pbl,ngrid
	74	integer nlay,nlev
	75
	76	logical first
	77	integer ipas
	78	save first,ipas
	79	!cFH/IM data first,ipas/.true.,0/
	80	data first,ipas/.false.,0/
	81	!$OMP THREADPRIVATE( first,ipas)
	82
	83	integer ig,k
	84
	85
	86	real ri,zrif,zalpha,zsm,zsn
	87	real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
	88
	89	real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
	90	real dtetadz(klon,klev+1)
	91	real m2cstat,mcstat,kmcstat
	92	real l(klon,klev+1)
	93	real,allocatable,save :: l0(:)
	94	!$OMP THREADPRIVATE(l0)
	95	real sq(klon),sqz(klon),zz(klon,klev+1)
	96	integer iter
	97
	98	real ric,rifc,b1,kap
	99	save ric,rifc,b1,kap
	100	data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
	101	!$OMP THREADPRIVATE(ric,rifc,b1,kap)
	102	real frif,falpha,fsm
	103	real fl,zzz,zl0,zq2,zn2
	104
	105	real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev) &
	106	& ,lyam(klon,klev),knyam(klon,klev) &
	107	& ,w2yam(klon,klev),t2yam(klon,klev)
	108	logical,save :: firstcall=.true.
	109	!$OMP THREADPRIVATE(firstcall)
	110	frif(ri)=0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
	111	falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
	112	fsm(ri)=1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
	113	fl(zzz,zl0,zq2,zn2)= &
	114	& max(min(l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig)) &
	115	& ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
	116
	117
	118	nlay=klev
	119	nlev=klev+1
	120
	121	if (firstcall) then
	122	allocate(l0(klon))
	123	firstcall=.false.
	124	endif
	125
	126
	127	if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.12)) then
	128	stop'probleme de coherence dans appel a MY'
	129	endif
	130
	131	ipas=ipas+1
	132
	133
	134	!c.......................................................................
	135	!c les increments verticaux
	136	!c.......................................................................
	137	!c
	138	!c!!!!! allerte !!!!!c
	139	!c!!!!! zlev n'est pas declare a nlev !!!!!c
	140	!c!!!!! ---->
	141	DO ig=1,ngrid
	142	zlev(ig,nlev)=zlay(ig,nlay) &
	143	& +( zlay(ig,nlay) - zlev(ig,nlev-1) )
	144	ENDDO
	145	!c!!!!! <----
	146	!c!!!!! allerte !!!!!c
	147	!c
	148	DO k=1,nlay
	149	DO ig=1,ngrid
	150	unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
	151	ENDDO
	152	ENDDO
	153	DO ig=1,ngrid
	154	unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
	155	ENDDO
	156	DO k=2,nlay
	157	DO ig=1,ngrid
	158	unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
	159	ENDDO
	160	ENDDO
	161	DO ig=1,ngrid
	162	unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
	163	ENDDO
	164	!c
	165	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	166	! Computing M^2, N^2, Richardson numbers, stability functions
	167	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	168
	169	do k=2,klev
	170	do ig=1,ngrid
	171	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
	172	m2(ig,k)=((u(ig,k)-u(ig,k-1))2+(v(ig,k)-v(ig,k-1))2) &
	173	& /(dz(ig,k)*dz(ig,k))
	174	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
	175	n2(ig,k)=g2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
	176	!c n2(ig,k)=0.
	177	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
	178	if (ri.lt.ric) then
	179	rif(ig,k)=frif(ri)
	180	else
	181	rif(ig,k)=rifc
	182	endif
	183	if(rif(ig,k).lt.0.16) then
	184	alpha(ig,k)=falpha(rif(ig,k))
	185	sm(ig,k)=fsm(rif(ig,k))
	186	else
	187	alpha(ig,k)=1.12
	188	sm(ig,k)=0.085
	189	endif
	190	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
	191	enddo
	192	enddo
	193
	194
	195	!c====================================================================
	196	!c Computing the mixing length
	197	!c====================================================================
	198
	199	!c Mise a jour de l0
	200	if (iflag_pbl==8.or.iflag_pbl==10) then
	201
	202	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	203	! Iterative computation of l0
	204	! This version is kept for iflag_pbl only for convergence
	205	! with NPv3.1 Cmip5 simulations
	206	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	207
	208	do ig=1,ngrid
	209	sq(ig)=1.e-10
	210	sqz(ig)=1.e-10
	211	enddo
	212	do k=2,klev-1
	213	do ig=1,ngrid
	214	zq=sqrt(q2(ig,k))
	215	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	216	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	217	enddo
	218	enddo
	219	do ig=1,ngrid
	220	l0(ig)=0.2*sqz(ig)/sq(ig)
	221	enddo
	222	do k=2,klev
	223	do ig=1,ngrid
	224	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	225	enddo
	226	enddo
	227	! print*,'L0 cas 8 ou 10 ',l0
	228
	229	else
	230
	231	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	232	! In all other case, the assymptotic mixing length l0 is imposed (100m)
	233	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	234
	235	l0(:)=150.
	236	do k=2,klev
	237	do ig=1,ngrid
	238	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	239	enddo
	240	enddo
	241	! print*,'L0 cas autres ',l0
	242
	243	endif
	244
	245
	246	!c====================================================================
	247	!c Yamada 2.0
	248	!c====================================================================
	249	if (iflag_pbl.eq.6) then
	250
	251	do k=2,klev
	252	q2(:,k)=l(:,k)*2zz(:,k)
	253	enddo
	254
	255
	256	else if (iflag_pbl.eq.7) then
	257	!c====================================================================
	258	!c Yamada 2.Fournier
	259	!c====================================================================
	260
	261	!c Calcul de l, km, au pas precedent
	262	do k=2,klev
	263	do ig=1,ngrid
	264	!c print*,'SMML=',sm(ig,k),l(ig,k)
	265	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	266	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	267	mpre(ig,k)=sqrt(m2(ig,k))
	268	!c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	269	enddo
	270	enddo
	271
	272	do k=2,klev-1
	273	do ig=1,ngrid
	274	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
	275	mcstat=sqrt(m2cstat)
	276
	277	!c print*,'M2 L=',k,mpre(ig,k),mcstat
	278	!c
	279	!c -----{puis on ecrit la valeur de q qui annule l'equation de m
	280	!c supposee en q3}
	281	!c
	282	IF (k.eq.2) THEN
	283	kmcstat=1.E+0 / mcstat &
	284	& ( unsdz(ig,k)kmpre(ig,k+1) &
	285	& *mpre(ig,k+1) &
	286	& +unsdz(ig,k-1) &
	287	& *cd(ig) &
	288	& ( sqrt(u(ig,3)2+v(ig,3)*2) &
	289	& -mcstat/unsdzdec(ig,k) &
	290	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2) &
	291	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	292	ELSE
	293	kmcstat=1.E+0 / mcstat &
	294	& ( unsdz(ig,k)kmpre(ig,k+1) &
	295	& *mpre(ig,k+1) &
	296	& +unsdz(ig,k-1)*kmpre(ig,k-1) &
	297	& *mpre(ig,k-1) ) &
	298	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	299	ENDIF
	300	!c print*,'T2 L=',k,tmp2
	301	tmp2=kmcstat &
	302	& /( sm(ig,k)/q2(ig,k) ) &
	303	& /l(ig,k)
	304	q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
	305	!c print*,'Q2 L=',k,q2(ig,k)
	306	!c
	307	enddo
	308	enddo
	309
	310	else if (iflag_pbl==8.or.iflag_pbl==9) then
	311	!c====================================================================
	312	!c Yamada 2.5 a la Didi
	313	!c====================================================================
	314
	315
	316	!c Calcul de l, km, au pas precedent
	317	do k=2,klev
	318	do ig=1,ngrid
	319	!c print*,'SMML=',sm(ig,k),l(ig,k)
	320	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	321	if (delta(ig,k).lt.1.e-20) then
	322	!c print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
	323	delta(ig,k)=1.e-20
	324	endif
	325	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	326	aa0= &
	327	& (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
	328	aa1= &
	329	& (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
	330	!c abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
	331	aa(ig,k)=aa1dt/(delta(ig,k)l(ig,k))
	332	!c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	333	qpre=sqrt(q2(ig,k))
	334	! if (iflag_pbl.eq.8 ) then
	335	if (aa(ig,k).gt.0.) then
	336	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
	337	else
	338	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	339	endif
	340	! else ! iflag_pbl=9
	341	! if (aa(ig,k)*qpre.gt.0.9) then
	342	! q2(ig,k)=(qpre10.)*2
	343	! else
	344	! q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	345	! endif
	346	! endif
	347	q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
	348	!c print,'Q2 L=',k,q2(ig,k),qpreqpre
	349	enddo
	350	enddo
	351
	352	else if (iflag_pbl>=10) then
	353
	354	! print*,'Schema mixte D'
	355	! print*,'Longueur ',l(:,:)
	356	do k=2,klev-1
	357	l(:,k)=max(l(:,k),1.)
	358	km(:,k)=l(:,k)sqrt(q2(:,k))sm(:,k)
	359	q2(:,k)=q2(:,k)+dtkm(:,k)m2(:,k)*(1.-rif(:,k))
	360	q2(:,k)=min(max(q2(:,k),1.e-10),1.e4)
	361	q2(:,k)=1./(1./sqrt(q2(:,k))+dt/(2l(:,k)b1))
	362	q2(:,k)=q2(:,k)*q2(:,k)
	363	enddo
	364
	365
	366	else
	367	stop'Cas nom prevu dans yamada4'
	368
	369	endif ! Fin du cas 8
	370
	371	!c print*,'OK8'
	372
	373	!c====================================================================
	374	!c Calcul des coefficients de mï¿œange
	375	!c====================================================================
	376	do k=2,klev
	377	!c print*,'k=',k
	378	do ig=1,ngrid
	379	!cabde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
	380	zq=sqrt(q2(ig,k))
	381	km(ig,k)=l(ig,k)zqsm(ig,k)
	382	kn(ig,k)=km(ig,k)*alpha(ig,k)
	383	kq(ig,k)=l(ig,k)zq0.2
	384	!c print*,'KML=',km(ig,k),kn(ig,k)
	385	enddo
	386	enddo
	387
	388	! Transport diffusif vertical de la TKE.
	389	if (iflag_pbl.ge.12) then
	390	! print*,'YAMADA VDIF'
	391	q2(:,1)=q2(:,2)
	392	call vdif_q2(dt,g,rconst,ngrid,plev,temp,kq,q2)
	393	endif
	394
	395	!c Traitement des cas noctrunes avec l'introduction d'une longueur
	396	!c minilale.
	397
	398	!c====================================================================
	399	!c Traitement particulier pour les cas tres stables.
	400	!c D'apres Holtslag Boville.
	401
	402	if (prt_level>1) THEN
	403	print*,'YAMADA4 0'
	404	endif !(prt_level>1) THEN
	405	do ig=1,ngrid
	406	coriol(ig)=1.e-4
	407	pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
	408	enddo
	409
	410	! print*,'pblhmin ',pblhmin
	411	!CTest a remettre 21 11 02
	412	!c test abd 13 05 02 if(0.eq.1) then
	413	if(1==1) then
	414	if(iflag_pbl==8.or.iflag_pbl==10) then
	415
	416	do k=2,klev
	417	do ig=1,ngrid
	418	if (teta(ig,2).gt.teta(ig,1)) then
	419	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	420	kmin=kapzlev(ig,k)qmin
	421	else
	422	kmin=-1. ! kmin n'est utilise que pour les SL stables.
	423	endif
	424	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
	425	!c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
	426	!c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
	427	kn(ig,k)=kmin
	428	km(ig,k)=kmin
	429	kq(ig,k)=kmin
	430	!c la longueur de melange est suposee etre l= kap z
	431	!c K=l q Sm d'ou q2=(K/l Sm)**2
	432	q2(ig,k)=(qmin/sm(ig,k))**2
	433	endif
	434	enddo
	435	enddo
	436
	437	else
	438
	439	do k=2,klev
	440	do ig=1,ngrid
	441	if (teta(ig,2).gt.teta(ig,1)) then
	442	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	443	kmin=kapzlev(ig,k)qmin
	444	else
	445	kmin=-1. ! kmin n'est utilise que pour les SL stables.
	446	endif
	447	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
	448	!c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
	449	!c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
	450	kn(ig,k)=kmin
	451	km(ig,k)=kmin
	452	kq(ig,k)=kmin
	453	!c la longueur de melange est suposee etre l= kap z
	454	!c K=l q Sm d'ou q2=(K/l Sm)**2
	455	sm(ig,k)=1.
	456	alpha(ig,k)=1.
	457	q2(ig,k)=min((qmin/sm(ig,k))**2,10.)
	458	zq=sqrt(q2(ig,k))
	459	km(ig,k)=l(ig,k)zqsm(ig,k)
	460	kn(ig,k)=km(ig,k)*alpha(ig,k)
	461	kq(ig,k)=l(ig,k)zq0.2
	462	endif
	463	enddo
	464	enddo
	465	endif
	466
	467	endif
	468
	469	if (prt_level>1) THEN
	470	print*,'YAMADA4 1'
	471	endif !(prt_level>1) THEN
	472	!c Diagnostique pour stokage
	473
	474	if(1.eq.0)then
	475	rino=rif
	476	smyam(1:ngrid,1)=0.
	477	styam(1:ngrid,1)=0.
	478	lyam(1:ngrid,1)=0.
	479	knyam(1:ngrid,1)=0.
	480	w2yam(1:ngrid,1)=0.
	481	t2yam(1:ngrid,1)=0.
	482
	483	smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
	484	styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
	485	lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
	486	knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
	487
	488	!c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
	489
	490	w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24 &
	491	& +lyam(1:ngrid,2:klev)5.17kn(1:ngrid,2:klev) &
	492	& *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
	493
	494	t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev) &
	495	& dtetadz(1:ngrid,2:klev)*2 &
	496	& /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
	497	endif
	498
	499	!c print*,'OKFIN'
	500	first=.false.
	501	return
	502	end
	503	SUBROUTINE vdif_q2(timestep,gravity,rconst,ngrid,plev,temp, &
	504	& kmy,q2)
	505	use dimphy
	506	IMPLICIT NONE
	507	!c.......................................................................
	508	#include "dimensions.h"
	509	!cccc#include "dimphy.h"
	510	!c.......................................................................
	511	!c
	512	!c dt : pas de temps
	513
	514	real plev(klon,klev+1)
	515	real temp(klon,klev)
	516	real timestep
	517	real gravity,rconst
	518	real kstar(klon,klev+1),zz
	519	real kmy(klon,klev+1)
	520	real q2(klon,klev+1)
	521	real deltap(klon,klev+1)
	522	real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
	523	integer ngrid
	524
	525	integer i,k
	526
	527	! print*,'RD=',rconst
	528	do k=1,klev
	529	do i=1,ngrid
	530	!c test
	531	! print*,'i,k',i,k
	532	! print*,'temp(i,k)=',temp(i,k)
	533	! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
	534	zz=(plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
	535	kstar(i,k)=0.125(kmy(i,k+1)+kmy(i,k))zz*zz &
	536	& /(plev(i,k)-plev(i,k+1))*timestep
	537	enddo
	538	enddo
	539
	540	do k=2,klev
	541	do i=1,ngrid
	542	deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
	543	enddo
	544	enddo
	545	do i=1,ngrid
	546	deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
	547	deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
	548	denom(i,klev+1)=deltap(i,klev+1)+kstar(i,klev)
	549	alpha(i,klev+1)=deltap(i,klev+1)*q2(i,klev+1)/denom(i,klev+1)
	550	beta(i,klev+1)=kstar(i,klev)/denom(i,klev+1)
	551	enddo
	552
	553	do k=klev,2,-1
	554	do i=1,ngrid
	555	denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))* &
	556	& kstar(i,k)+kstar(i,k-1)
	557	!c correction d'un bug 10 01 2001
	558	alpha(i,k)=(q2(i,k)*deltap(i,k) &
	559	& +kstar(i,k)*alpha(i,k+1))/denom(i,k)
	560	beta(i,k)=kstar(i,k-1)/denom(i,k)
	561	enddo
	562	enddo
	563
	564	!c Si on recalcule q2(1)
	565	if(1.eq.0) then
	566	do i=1,ngrid
	567	denom(i,1)=deltap(i,1)+(1-beta(i,2))*kstar(i,1)
	568	q2(i,1)=(q2(i,1)*deltap(i,1) &
	569	& +kstar(i,1)*alpha(i,2))/denom(i,1)
	570	enddo
	571	endif
	572	!c sinon, on peut sauter cette boucle...
	573
	574	do k=2,klev+1
	575	do i=1,ngrid
	576	q2(i,k)=alpha(i,k)+beta(i,k)*q2(i,k-1)
	577	enddo
	578	enddo
	579
	580	return
	581	end
	582	SUBROUTINE vdif_q2e(timestep,gravity,rconst,ngrid, &
	583	& plev,temp,kmy,q2)
	584	use dimphy
	585	IMPLICIT NONE
	586	!c.......................................................................
	587	#include "dimensions.h"
	588	!cccc#include "dimphy.h"
	589	!c.......................................................................
	590	!c
	591	!c dt : pas de temps
	592
	593	real plev(klon,klev+1)
	594	real temp(klon,klev)
	595	real timestep
	596	real gravity,rconst
	597	real kstar(klon,klev+1),zz
	598	real kmy(klon,klev+1)
	599	real q2(klon,klev+1)
	600	real deltap(klon,klev+1)
	601	real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
	602	integer ngrid
	603
	604	integer i,k
	605
	606	do k=1,klev
	607	do i=1,ngrid
	608	zz=(plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
	609	kstar(i,k)=0.125(kmy(i,k+1)+kmy(i,k))zz*zz &
	610	& /(plev(i,k)-plev(i,k+1))*timestep
	611	enddo
	612	enddo
	613
	614	do k=2,klev
	615	do i=1,ngrid
	616	deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
	617	enddo
	618	enddo
	619	do i=1,ngrid
	620	deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
	621	deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
	622	enddo
	623
	624	do k=klev,2,-1
	625	do i=1,ngrid
	626	q2(i,k)=q2(i,k)+ &
	627	& ( kstar(i,k)*(q2(i,k+1)-q2(i,k)) &
	628	& -kstar(i,k-1)*(q2(i,k)-q2(i,k-1)) ) &
	629	& /deltap(i,k)
	630	enddo
	631	enddo
	632
	633	do i=1,ngrid
	634	q2(i,1)=q2(i,1)+ &
	635	& ( kstar(i,1)*(q2(i,2)-q2(i,1)) &
	636	& ) &
	637	& /deltap(i,1)
	638	q2(i,klev+1)=q2(i,klev+1)+ &
	639	& ( &
	640	& -kstar(i,klev)*(q2(i,klev+1)-q2(i,klev)) ) &
	641	& /deltap(i,klev+1)
	642	enddo
	643
	644	return
	645	end

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

Download in other formats:

Original Format