Context Navigation

yamada4.F @ 1461

Last change on this file since 1461 was 814, checked in by slebonnois, 12 years ago
SL: corrections de bugs pour Titan, + suppression de la dependance circulaire dans mod_grid_phy_lmdz.F90, + petite inversion sur la ligne 101 de create_make_gcm pour que disvert compile du premier coup.
File size: 15.8 KB

Rev	Line
[3]	1	!
	2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/yamada4.F,v 1.1 2004/06/22 11:45:36 lmdzadmin Exp $
	3	!
	4	SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
	5	c s ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
	6	s ,zlev,zlay,u,v,teta,cd,km,kn,kq,ustar
	7	s ,iflag_pbl)
[102]	8	c.......................................................................
	9	use dimphy
[3]	10	IMPLICIT NONE
	11	#include "dimensions.h"
[814]	12	#include "YOMCST.h"
[3]	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 : MY 2.5
	40	c iflag_pbl=9 : un test ?
	41
	42	c.......................................................................
	43	REAL dt,g,rconst
	44	real plev(klon,klev+1),temp(klon,klev)
	45	real ustar(klon)
	46	real kmin,qmin,pblhmin(klon),coriol(klon)
	47	REAL zlev(klon,klev+1)
	48	REAL zlay(klon,klev)
	49	REAL u(klon,klev)
	50	REAL v(klon,klev)
	51	REAL teta(klon,klev)
	52	REAL cd(klon)
[102]	53	REAL qpre
[3]	54	REAL unsdz(klon,klev)
	55	REAL unsdzdec(klon,klev+1)
	56
	57	REAL km(klon,klev+1)
	58	REAL kmpre(klon,klev+1),tmp2
	59	REAL mpre(klon,klev+1)
	60	REAL kn(klon,klev+1)
	61	REAL kq(klon,klev+1)
	62	real ff(klon,klev+1),delta(klon,klev+1)
	63	real aa(klon,klev+1),aa0,aa1
	64	integer iflag_pbl,ngrid
	65
	66
	67	integer nlay,nlev
	68
	69	logical first
	70	integer ipas
	71	save first,ipas
	72	data first,ipas/.true.,0/
	73
	74
	75	integer ig,k
	76
	77
	78	real ri,zrif,zalpha,zsm,zsn
	79	real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
	80
	81	real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
	82	real dtetadz(klon,klev+1)
	83	real m2cstat,mcstat,kmcstat
[102]	84	real l(klon,klev+1)
	85	real,save,allocatable :: l0(:)
[3]	86	c ATTENTION! mis ici car j'ai enlevé q2 des arguments...
	87	c sinon, c'est au-dessus que ça se passe...
[102]	88	REAL,save,allocatable :: q2(:,:)
[3]	89
	90	real sq(klon),sqz(klon),zz(klon,klev+1)
	91	integer iter
	92
	93	real ric,rifc,b1,kap
	94	save ric,rifc,b1,kap
	95	data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
	96
	97	real frif,falpha,fsm
	98	real fl,zzz,zl0,zq2,zn2
	99
[102]	100	c real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
	101	c s ,lyam(klon,klev),knyam(klon,klev)
	102	c s ,w2yam(klon,klev),t2yam(klon,klev)
[3]	103
	104	frif(ri)=0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
	105	falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
	106	fsm(ri)=1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
	107	fl(zzz,zl0,zq2,zn2)=
	108	s max(min(l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
	109	s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
	110
	111	if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
	112	stop'probleme de coherence dans appel a MY'
	113	endif
	114
	115	c===================================
[102]	116	c INITIALISATIONS
	117	nlay=klev
	118	nlev=klev+1
	119
[3]	120	if (first) then
[102]	121	allocate(l0(klon))
	122	allocate(q2(klon,klevp1))
	123
	124	c (surtout pour k=1, à cause diagnostiques...)
[3]	125	q2 = 0.
	126	endif
	127	c===================================
	128
	129	ipas=ipas+1
	130	if (0.eq.1.and.first) then
	131	do ig=1,1000
	132	ri=(ig-800.)/500.
	133	if (ri.lt.ric) then
	134	zrif=frif(ri)
	135	else
	136	zrif=rifc
	137	endif
	138	if(zrif.lt.0.16) then
	139	zalpha=falpha(zrif)
	140	zsm=fsm(zrif)
	141	else
	142	zalpha=1.12
	143	zsm=0.085
	144	endif
	145	c print*,ri,rif,zalpha,zsm
	146	enddo
	147	endif
	148
	149	c.......................................................................
	150	c les increments verticaux
	151	c.......................................................................
	152	c
	153	c!!!!! allerte !!!!!c
	154	c!!!!! zlev n'est pas declare a nlev !!!!!c
	155	c!!!!! ---->
	156	DO ig=1,ngrid
	157	zlev(ig,nlev)=zlay(ig,nlay)
	158	& +( zlay(ig,nlay) - zlev(ig,nlev-1) )
	159	ENDDO
	160	c!!!!! <----
	161	c!!!!! allerte !!!!!c
	162	c
	163	DO k=1,nlay
	164	DO ig=1,ngrid
	165	unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
	166	ENDDO
	167	ENDDO
	168	DO ig=1,ngrid
	169	unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
	170	ENDDO
	171	DO k=2,nlay
	172	DO ig=1,ngrid
	173	unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
	174	ENDDO
	175	ENDDO
	176	DO ig=1,ngrid
	177	unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
	178	ENDDO
	179	c
	180	c.......................................................................
	181
	182	c===================================
	183	c INITIALISATIONS (surtout pour k=1, à cause diagnostiques...)
	184	dz = 0.
	185	m2 = 0.
	186	dtetadz = 0.
	187	n2 = 0.
	188	rif = 0.
	189	alpha = 0.
	190	sm = 0.
	191	zz = 0.
	192	l = 0.
	193	c===================================
	194	do k=2,klev
	195	do ig=1,ngrid
	196	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
	197	m2(ig,k)=((u(ig,k)-u(ig,k-1))2+(v(ig,k)-v(ig,k-1))2)
	198	s /(dz(ig,k)*dz(ig,k))
	199	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
	200	n2(ig,k)=g2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
	201	c n2(ig,k)=0.
	202	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
	203	if (ri.lt.ric) then
	204	rif(ig,k)=frif(ri)
	205	else
	206	rif(ig,k)=rifc
	207	endif
	208	if(rif(ig,k).lt.0.16) then
	209	alpha(ig,k)=falpha(rif(ig,k))
	210	sm(ig,k)=fsm(rif(ig,k))
	211	else
	212	alpha(ig,k)=1.12
	213	sm(ig,k)=0.085
	214	endif
	215	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
	216	c print,'RIF L=',k,rif(ig,k),rialpha(ig,k)
	217
	218
	219	enddo
	220	enddo
	221
	222
	223	c====================================================================
	224	c Au premier appel, on determine l et q2 de facon iterative.
	225	c iterration pour determiner la longueur de melange
	226
	227
	228	if (first.or.iflag_pbl.eq.6) then
	229	do ig=1,ngrid
	230	l0(ig)=10.
	231	enddo
	232	do k=2,klev-1
	233	do ig=1,ngrid
	234	l(ig,k)=l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
	235	enddo
	236	enddo
	237
	238	do iter=1,10
	239	do ig=1,ngrid
	240	sq(ig)=1.e-10
	241	sqz(ig)=1.e-10
	242	enddo
	243	do k=2,klev-1
	244	do ig=1,ngrid
	245	q2(ig,k)=l(ig,k)*2zz(ig,k)
	246	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	247	zq=sqrt(q2(ig,k))
	248	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	249	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	250	enddo
	251	enddo
	252	do ig=1,ngrid
	253	l0(ig)=0.2*sqz(ig)/sq(ig)
	254	c l0(ig)=30.
	255	enddo
	256	c print*,'ITER=',iter,' L0=',l0
	257
	258	enddo
	259
	260	c print*,'Fin de l initialisation de q2 et l0'
	261
	262	endif ! first
	263
	264	c====================================================================
	265	c Calcul de la longueur de melange.
	266	c====================================================================
	267
	268	c Mise a jour de l0
	269	do ig=1,ngrid
	270	sq(ig)=1.e-10
	271	sqz(ig)=1.e-10
	272	enddo
	273	do k=2,klev-1
	274	do ig=1,ngrid
	275	zq=sqrt(q2(ig,k))
	276	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	277	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	278	enddo
	279	enddo
	280	do ig=1,ngrid
	281	l0(ig)=0.2*sqz(ig)/sq(ig)
	282	c l0(ig)=30.
	283	enddo
	284	c print*,'ITER=',iter,' L0=',l0
	285	c calcul de l(z)
	286	do k=2,klev
	287	do ig=1,ngrid
	288	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	289	if(first) then
	290	q2(ig,k)=l(ig,k)*2zz(ig,k)
	291	endif
	292	enddo
	293	enddo
	294
	295	c====================================================================
	296	c Yamada 2.0
	297	c====================================================================
	298	if (iflag_pbl.eq.6) then
	299
	300	do k=2,klev
	301	do ig=1,ngrid
	302	q2(ig,k)=l(ig,k)*2zz(ig,k)
	303	enddo
	304	enddo
	305
	306
	307	else if (iflag_pbl.eq.7) then
	308	c====================================================================
	309	c Yamada 2.Fournier
	310	c====================================================================
	311
	312	c Calcul de l, km, au pas precedent
	313	do k=2,klev
	314	do ig=1,ngrid
	315	c print*,'SMML=',sm(ig,k),l(ig,k)
	316	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	317	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	318	mpre(ig,k)=sqrt(m2(ig,k))
	319	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	320	enddo
	321	enddo
	322
	323	do k=2,klev-1
	324	do ig=1,ngrid
	325	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
	326	mcstat=sqrt(m2cstat)
	327
	328	c print*,'M2 L=',k,mpre(ig,k),mcstat
	329	c
	330	c -----{puis on ecrit la valeur de q qui annule l'equation de m
	331	c supposee en q3}
	332	c
	333	IF (k.eq.2) THEN
	334	kmcstat=1.E+0 / mcstat
	335	& ( unsdz(ig,k)kmpre(ig,k+1)
	336	& *mpre(ig,k+1)
	337	& +unsdz(ig,k-1)
	338	& *cd(ig)
	339	& ( sqrt(u(ig,3)2+v(ig,3)*2)
	340	& -mcstat/unsdzdec(ig,k)
	341	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
	342	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	343	ELSE
	344	kmcstat=1.E+0 / mcstat
	345	& ( unsdz(ig,k)kmpre(ig,k+1)
	346	& *mpre(ig,k+1)
	347	& +unsdz(ig,k-1)*kmpre(ig,k-1)
	348	& *mpre(ig,k-1) )
	349	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	350	ENDIF
	351	c print*,'T2 L=',k,tmp2
	352	tmp2=kmcstat
	353	& /( sm(ig,k)/q2(ig,k) )
	354	& /l(ig,k)
	355	q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
	356	c print*,'Q2 L=',k,q2(ig,k)
	357	c
	358	enddo
	359	enddo
	360
	361	else if (iflag_pbl.ge.8) then
	362	c====================================================================
	363	c Yamada 2.5 a la Didi
	364	c====================================================================
	365
	366
	367	c Calcul de l, km, au pas precedent
	368	do k=2,klev
	369	do ig=1,ngrid
	370	c print*,'SMML=',sm(ig,k),l(ig,k)
	371	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	372	if (delta(ig,k).lt.1.e-20) then
	373	c print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
	374	delta(ig,k)=1.e-20
	375	endif
	376	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	377	aa0=
	378	s (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
	379	aa1=
	380	s (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
	381	c abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
	382	aa(ig,k)=aa1dt/(delta(ig,k)l(ig,k))
	383	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	384	qpre=sqrt(q2(ig,k))
	385	if (iflag_pbl.eq.8 ) then
	386	if (aa(ig,k).gt.0.) then
	387	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
	388	else
	389	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	390	endif
	391	else ! iflag_pbl=9
	392	if (aa(ig,k)*qpre.gt.0.9) then
	393	q2(ig,k)=(qpre10.)*2
	394	else
	395	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	396	endif
	397	endif
	398	q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
	399	c print,'Q2 L=',k,q2(ig,k),qpreqpre
	400	enddo
	401	enddo
	402
	403	endif ! Fin du cas 8
	404
	405	c print*,'OK8'
	406
	407	c====================================================================
	408	c Calcul des coefficients de mélange
	409	c====================================================================
	410	do k=2,klev
	411	c print*,'k=',k
	412	do ig=1,ngrid
	413	cabde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
	414	zq=sqrt(q2(ig,k))
	415	km(ig,k)=l(ig,k)zqsm(ig,k)
	416	kn(ig,k)=km(ig,k)*alpha(ig,k)
	417	kq(ig,k)=l(ig,k)zq0.2
	418	c print*,'KML=',km(ig,k),kn(ig,k)
	419	enddo
	420	enddo
	421
	422	c if (iflag_pbl.ge.7..and.0.eq.1) then
	423	c q2(:,1)=q2(:,2)
	424	c call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
	425	c endif
	426
	427	c Traitement des cas noctrunes avec l'introduction d'une longueur
	428	c minilale.
	429
	430	c====================================================================
	431	c Traitement particulier pour les cas tres stables.
	432	c D'apres Holtslag Boville.
	433
	434	c print*,'YAMADA4 0'
	435
	436	do ig=1,ngrid
[808]	437	coriol(ig)=1.e-4*86400/RDAY !! scaling... should be checked
[814]	438	pblhmin(ig)=0.07*ustar(ig)/
	439	. max(abs(coriol(ig)),2.546e-5*86400/RDAY)
[3]	440	enddo
[808]	441	c if (first) then
	442	c print*,'A REVOIR!! coriol ?? pblhmin ',pblhmin
	443	c endif
[3]	444	CTest a remettre 21 11 02
	445	c test abd 13 05 02 if(0.eq.1) then
	446	if(1.eq.1) then
	447	do k=2,klev
	448	do ig=1,klon
	449	if (teta(ig,2).gt.teta(ig,1)) then
	450	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	451	kmin=kapzlev(ig,k)qmin
	452	else
	453	kmin=-1. ! kmin n'est utilise que pour les SL stables.
	454	endif
	455	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
	456	c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
	457	c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
	458	kn(ig,k)=kmin
	459	km(ig,k)=kmin
	460	kq(ig,k)=kmin
	461	c la longueur de melange est suposee etre l= kap z
	462	c K=l q Sm d'ou q2=(K/l Sm)**2
	463	q2(ig,k)=(qmin/sm(ig,k))**2
	464	endif
	465	enddo
	466	enddo
	467	endif
	468
	469	c print*,'YAMADA4 1'
	470
	471	c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
	472
[102]	473	c if(1.eq.0)then
	474	c w2yam=q2(:,1:klev)*0.24
	475	c s +lyam(:,1:klev)5.17kn(:,1:klev)*n2(:,1:klev)
	476	c s /sqrt(q2(:,1:klev))
	477	c
	478	c t2yam=9.1kn(:,1:klev)dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))
	479	c s *lyam(:,1:klev)
	480	c endif
[3]	481
	482	c print*,'OKFIN'
	483	first=.false.
	484	return
	485	end

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

Download in other formats:

Original Format