Context Navigation

yamada4.F @ 3556

Last change on this file since 3556 was 2535, checked in by emillour, 4 years ago
Venus GCM: Update yamada4 to use q2 computed from previous time step (or read from startphy) so that it inow conforms to 1+1=2. Some code tidying in clmain and yamada4 along the way. VB + EM
File size: 15.3 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
[2535]	5	s ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
[3]	6	s ,iflag_pbl)
[101]	7	c.......................................................................
[2535]	8	use dimphy, only: klon, klev
	9	use turb_mod, only: l0
	10
	11	#ifdef CPP_XIOS
	12	use xios_output_mod, only: send_xios_field
	13	#endif
	14
[3]	15	IMPLICIT NONE
	16	c.......................................................................
	17	c
	18	c dt : pas de temps
	19	c g : g
	20	c zlev : altitude a chaque niveau (interface inferieure de la couche
	21	c de meme indice)
	22	c zlay : altitude au centre de chaque couche
	23	c u,v : vitesse au centre de chaque couche
	24	c (en entree : la valeur au debut du pas de temps)
	25	c teta : temperature potentielle au centre de chaque couche
	26	c (en entree : la valeur au debut du pas de temps)
	27	c cd : cdrag
	28	c (en entree : la valeur au debut du pas de temps)
	29	c q2 : $q^2$ au bas de chaque couche
	30	c (en entree : la valeur au debut du pas de temps)
	31	c (en sortie : la valeur a la fin du pas de temps)
	32	c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
	33	c couche)
	34	c (en sortie : la valeur a la fin du pas de temps)
	35	c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
	36	c (en sortie : la valeur a la fin du pas de temps)
	37	c
	38	c iflag_pbl doit valoir entre 6 et 9
	39	c l=6, on prend systematiquement une longueur d'equilibre
	40	c iflag_pbl=6 : MY 2.0
	41	c iflag_pbl=7 : MY 2.0.Fournier
	42	c iflag_pbl=8 : MY 2.5
	43	c iflag_pbl=9 : un test ?
	44
	45	c.......................................................................
[2535]	46	REAL, INTENT(IN) :: dt,g,rconst
	47	REAL, INTENT(IN) :: plev(klon,klev+1),temp(klon,klev)
	48	REAL, INTENT(IN) :: ustar(klon)
	49	REAL, INTENT(INOUT) :: zlev(klon,klev+1)
	50	REAL, INTENT(IN) :: zlay(klon,klev)
	51	REAL, INTENT(IN) :: u(klon,klev)
	52	REAL, INTENT(IN) :: v(klon,klev)
	53	REAL, INTENT(IN) :: teta(klon,klev)
	54	REAL, INTENT(IN) :: cd(klon)
	55	INTEGER, INTENT(IN) :: iflag_pbl,ngrid
	56
	57	REAL, INTENT(OUT) :: km(klon,klev+1)
	58	REAL, INTENT(OUT) :: kn(klon,klev+1)
	59	REAL, INTENT(OUT) :: kq(klon,klev+1)
	60
	61	REAL, INTENT(INOUT) :: q2(klon,klev+1)
	62	c Variables locales:
[3]	63	real kmin,qmin,pblhmin(klon),coriol(klon)
[101]	64	REAL qpre
[3]	65	REAL unsdz(klon,klev)
	66	REAL unsdzdec(klon,klev+1)
	67
	68	REAL kmpre(klon,klev+1),tmp2
	69	REAL mpre(klon,klev+1)
[2535]	70
[3]	71	real ff(klon,klev+1),delta(klon,klev+1)
	72	real aa(klon,klev+1),aa0,aa1
[2535]	73
[3]	74	integer nlay,nlev
	75
[2535]	76	logical,save :: first=.false. ! not neede any more
	77	integer,save :: ipas=0
[3]	78
	79
	80	integer ig,k
	81
	82
	83	real ri,zrif,zalpha,zsm,zsn
	84	real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
	85
	86	real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
	87	real dtetadz(klon,klev+1)
	88	real m2cstat,mcstat,kmcstat
[101]	89	real l(klon,klev+1)
[3]	90
	91	real sq(klon),sqz(klon),zz(klon,klev+1)
	92	integer iter
	93
[2534]	94	real,save :: ric=0.195
	95	real,save :: rifc=0.191
	96	real,save :: b1=16.6
	97	real,save :: kap=0.4
[3]	98
	99	real frif,falpha,fsm
	100	real fl,zzz,zl0,zq2,zn2
[2535]	101
[101]	102	c real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
	103	c s ,lyam(klon,klev),knyam(klon,klev)
	104	c s ,w2yam(klon,klev),t2yam(klon,klev)
[3]	105
	106	frif(ri)=0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
	107	falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
	108	fsm(ri)=1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
	109	fl(zzz,zl0,zq2,zn2)=
	110	s max(min(l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
	111	s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
	112
[2535]	113
[3]	114	if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
	115	stop'probleme de coherence dans appel a MY'
	116	endif
	117
	118	c===================================
[101]	119	c INITIALISATIONS
	120	nlay=klev
	121	nlev=klev+1
	122
[2535]	123	c if (first) then
	124	c IF (.not.ALLOCATED(l0)) allocate(l0(klon))
	125	c IF (.not.ALLOCATED(q2)) allocate(q2(klon,klev+1))
	126	c endif
[3]	127	c===================================
[2535]	128
[3]	129	ipas=ipas+1
[2535]	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
[3]	145	c print*,ri,rif,zalpha,zsm
[2535]	146	c enddo
	147	c endif
	148
[3]	149	c.......................................................................
	150	c les increments verticaux
	151	c.......................................................................
	152	c
	153	DO k=1,nlay
[2535]	154	DO ig=1,ngrid
[3]	155	unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
[2535]	156	ENDDO
[3]	157	ENDDO
[2535]	158	DO ig=1,ngrid
[3]	159	unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
[2535]	160	ENDDO
[3]	161	DO k=2,nlay
[2535]	162	DO ig=1,ngrid
[3]	163	unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
[2535]	164	ENDDO
[3]	165	ENDDO
[2535]	166	DO ig=1,ngrid
[3]	167	unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
[2535]	168	ENDDO
[3]	169	c
	170	c.......................................................................
[2535]	171
[3]	172	c===================================
	173	c INITIALISATIONS (surtout pour k=1, à cause diagnostiques...)
[2535]	174	dz(:,:) = 0.
	175	m2(:,:) = 0.
	176	dtetadz(:,:) = 0.
	177	n2(:,:) = 0.
	178	rif(:,:) = 0.
	179	alpha(:,:) = 0.
	180	sm(:,:) = 0.
	181	zz(:,:) = 0.
	182	l(:,:) = 0.
	183	km(:,:) = 0.
	184	kn(:,:) = 0.
	185	kmpre(:,:) = 0.
	186	mpre(:,:) = 0.
	187
[3]	188	c===================================
	189	do k=2,klev
[2535]	190	do ig=1,ngrid
[3]	191	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
	192	m2(ig,k)=((u(ig,k)-u(ig,k-1))2+(v(ig,k)-v(ig,k-1))2)
	193	s /(dz(ig,k)*dz(ig,k))
	194	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
	195	n2(ig,k)=g2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
	196	c n2(ig,k)=0.
	197	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
	198	if (ri.lt.ric) then
	199	rif(ig,k)=frif(ri)
	200	else
	201	rif(ig,k)=rifc
	202	endif
	203	if(rif(ig,k).lt.0.16) then
	204	alpha(ig,k)=falpha(rif(ig,k))
	205	sm(ig,k)=fsm(rif(ig,k))
	206	else
	207	alpha(ig,k)=1.12
	208	sm(ig,k)=0.085
	209	endif
	210	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
	211	c print,'RIF L=',k,rif(ig,k),rialpha(ig,k)
	212
[2535]	213	enddo
[3]	214	enddo
	215
	216	c====================================================================
	217	c Au premier appel, on determine l et q2 de facon iterative.
	218	c iterration pour determiner la longueur de melange
	219
	220
	221	if (first.or.iflag_pbl.eq.6) then
	222	do ig=1,ngrid
	223	l0(ig)=10.
	224	enddo
	225	do k=2,klev-1
	226	do ig=1,ngrid
	227	l(ig,k)=l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
	228	enddo
	229	enddo
	230
	231	do iter=1,10
	232	do ig=1,ngrid
	233	sq(ig)=1.e-10
	234	sqz(ig)=1.e-10
	235	enddo
	236	do k=2,klev-1
	237	do ig=1,ngrid
	238	q2(ig,k)=l(ig,k)*2zz(ig,k)
	239	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	240	zq=sqrt(q2(ig,k))
	241	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	242	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	243	enddo
	244	enddo
	245	do ig=1,ngrid
	246	l0(ig)=0.2*sqz(ig)/sq(ig)
	247	c l0(ig)=30.
	248	enddo
	249	c print*,'ITER=',iter,' L0=',l0
	250
	251	enddo
	252
	253	c print*,'Fin de l initialisation de q2 et l0'
[2535]	254	endif ! of if (first.or.iflag_pbl.eq.6)
	255
[3]	256	c====================================================================
	257	c Calcul de la longueur de melange.
	258	c====================================================================
	259
	260	c Mise a jour de l0
	261	do ig=1,ngrid
	262	sq(ig)=1.e-10
	263	sqz(ig)=1.e-10
	264	enddo
	265	do k=2,klev-1
	266	do ig=1,ngrid
	267	zq=sqrt(q2(ig,k))
	268	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
	269	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
	270	enddo
	271	enddo
	272	do ig=1,ngrid
	273	l0(ig)=0.2*sqz(ig)/sq(ig)
	274	c l0(ig)=30.
	275	enddo
	276	c print*,'ITER=',iter,' L0=',l0
	277	c calcul de l(z)
	278	do k=2,klev
	279	do ig=1,ngrid
	280	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
	281	if(first) then
	282	q2(ig,k)=l(ig,k)*2zz(ig,k)
	283	endif
	284	enddo
	285	enddo
	286	c====================================================================
	287	c Yamada 2.0
	288	c====================================================================
	289	if (iflag_pbl.eq.6) then
	290
	291	do k=2,klev
	292	do ig=1,ngrid
	293	q2(ig,k)=l(ig,k)*2zz(ig,k)
	294	enddo
	295	enddo
	296
	297
	298	else if (iflag_pbl.eq.7) then
	299	c====================================================================
	300	c Yamada 2.Fournier
	301	c====================================================================
	302
	303	c Calcul de l, km, au pas precedent
	304	do k=2,klev
	305	do ig=1,ngrid
	306	c print*,'SMML=',sm(ig,k),l(ig,k)
	307	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	308	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	309	mpre(ig,k)=sqrt(m2(ig,k))
	310	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	311	enddo
	312	enddo
	313
	314	do k=2,klev-1
	315	do ig=1,ngrid
	316	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
	317	mcstat=sqrt(m2cstat)
	318
	319	c print*,'M2 L=',k,mpre(ig,k),mcstat
	320	c
	321	c -----{puis on ecrit la valeur de q qui annule l'equation de m
	322	c supposee en q3}
	323	c
	324	IF (k.eq.2) THEN
	325	kmcstat=1.E+0 / mcstat
	326	& ( unsdz(ig,k)kmpre(ig,k+1)
	327	& *mpre(ig,k+1)
	328	& +unsdz(ig,k-1)
	329	& *cd(ig)
	330	& ( sqrt(u(ig,3)2+v(ig,3)*2)
	331	& -mcstat/unsdzdec(ig,k)
	332	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
	333	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	334	ELSE
	335	kmcstat=1.E+0 / mcstat
	336	& ( unsdz(ig,k)kmpre(ig,k+1)
	337	& *mpre(ig,k+1)
	338	& +unsdz(ig,k-1)*kmpre(ig,k-1)
	339	& *mpre(ig,k-1) )
	340	& /( unsdz(ig,k)+unsdz(ig,k-1) )
	341	ENDIF
	342	c print*,'T2 L=',k,tmp2
	343	tmp2=kmcstat
	344	& /( sm(ig,k)/q2(ig,k) )
	345	& /l(ig,k)
	346	q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
	347	c print*,'Q2 L=',k,q2(ig,k)
	348	c
	349	enddo
	350	enddo
	351
	352	else if (iflag_pbl.ge.8) then
	353	c====================================================================
	354	c Yamada 2.5 a la Didi
	355	c====================================================================
	356
	357
	358	c Calcul de l, km, au pas precedent
	359	do k=2,klev
	360	do ig=1,ngrid
	361	c print*,'SMML=',sm(ig,k),l(ig,k)
	362	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
	363	if (delta(ig,k).lt.1.e-20) then
	364	c print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
	365	delta(ig,k)=1.e-20
	366	endif
	367	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
	368	aa0=
	369	s (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
	370	aa1=
	371	s (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
	372	c abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
	373	aa(ig,k)=aa1dt/(delta(ig,k)l(ig,k))
	374	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
	375	qpre=sqrt(q2(ig,k))
	376	if (iflag_pbl.eq.8 ) then
	377	if (aa(ig,k).gt.0.) then
	378	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
	379	else
	380	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	381	endif
	382	else ! iflag_pbl=9
	383	if (aa(ig,k)*qpre.gt.0.9) then
	384	q2(ig,k)=(qpre10.)*2
	385	else
	386	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
	387	endif
	388	endif
	389	q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
	390	c print,'Q2 L=',k,q2(ig,k),qpreqpre
	391	enddo
	392	enddo
	393
	394	endif ! Fin du cas 8
	395	c====================================================================
	396	c Calcul des coefficients de mélange
	397	c====================================================================
	398	do k=2,klev
	399	c print*,'k=',k
	400	do ig=1,ngrid
	401	cabde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
	402	zq=sqrt(q2(ig,k))
	403	km(ig,k)=l(ig,k)zqsm(ig,k)
	404	kn(ig,k)=km(ig,k)*alpha(ig,k)
	405	kq(ig,k)=l(ig,k)zq0.2
	406	c print*,'KML=',km(ig,k),kn(ig,k)
	407	enddo
	408	enddo
	409
	410	c if (iflag_pbl.ge.7..and.0.eq.1) then
	411	c q2(:,1)=q2(:,2)
	412	c call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
	413	c endif
	414
	415	c Traitement des cas noctrunes avec l'introduction d'une longueur
	416	c minilale.
	417
	418	c====================================================================
	419	c Traitement particulier pour les cas tres stables.
	420	c D'apres Holtslag Boville.
	421
	422	c print*,'YAMADA4 0'
	423
	424	do ig=1,ngrid
	425	coriol(ig)=1.e-4
	426	pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
	427	enddo
	428	if (first) then
	429	print*,'A REVOIR!! coriol ?? pblhmin ',pblhmin
	430	endif
	431	CTest a remettre 21 11 02
	432	c test abd 13 05 02 if(0.eq.1) then
	433	if(1.eq.1) then
	434	do k=2,klev
	435	do ig=1,klon
	436	if (teta(ig,2).gt.teta(ig,1)) then
	437	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
	438	kmin=kapzlev(ig,k)qmin
	439	else
	440	kmin=-1. ! kmin n'est utilise que pour les SL stables.
	441	endif
	442	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
	443	c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
	444	c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
	445	kn(ig,k)=kmin
	446	km(ig,k)=kmin
	447	kq(ig,k)=kmin
	448	c la longueur de melange est suposee etre l= kap z
	449	c K=l q Sm d'ou q2=(K/l Sm)**2
	450	q2(ig,k)=(qmin/sm(ig,k))**2
	451	endif
	452	enddo
	453	enddo
	454	endif
	455
	456	c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
	457
[101]	458	c if(1.eq.0)then
	459	c w2yam=q2(:,1:klev)*0.24
	460	c s +lyam(:,1:klev)5.17kn(:,1:klev)*n2(:,1:klev)
	461	c s /sqrt(q2(:,1:klev))
	462	c
	463	c t2yam=9.1kn(:,1:klev)dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))
	464	c s *lyam(:,1:klev)
	465	c endif
[3]	466
	467	c print*,'OKFIN'
	468	first=.false.
[2535]	469
[3]	470	end

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

Download in other formats:

Original Format