Context Navigation

vdif_kcay.F @ 3556

Last change on this file since 3556 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

File size: 25.1 KB

Rev	Line
[3]	1	!
	2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/vdif_kcay.F,v 1.1 2004/06/22 11:45:36 lmdzadmin Exp $
	3	!
	4	SUBROUTINE vdif_kcay(ngrid,dt,g,rconst,plev,temp
	5	s ,zlev,zlay,u,v,teta,cd,q2,q2diag,km,kn,ustar
	6	s ,l_mix)
[101]	7	c.......................................................................
	8	use dimphy
[3]	9	IMPLICIT NONE
	10	c.......................................................................
	11	c
	12	c dt : pas de temps
	13	c g : g
	14	c zlev : altitude a chaque niveau (interface inferieure de la couche
	15	c de meme indice)
	16	c zlay : altitude au centre de chaque couche
	17	c u,v : vitesse au centre de chaque couche
	18	c (en entree : la valeur au debut du pas de temps)
	19	c teta : temperature potentielle au centre de chaque couche
	20	c (en entree : la valeur au debut du pas de temps)
	21	c cd : cdrag
	22	c (en entree : la valeur au debut du pas de temps)
	23	c q2 : $q^2$ au bas de chaque couche
	24	c (en entree : la valeur au debut du pas de temps)
	25	c (en sortie : la valeur a la fin du pas de temps)
	26	c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
	27	c couche)
	28	c (en sortie : la valeur a la fin du pas de temps)
	29	c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
	30	c (en sortie : la valeur a la fin du pas de temps)
	31	c
	32	c.......................................................................
	33	REAL dt,g,rconst
	34	real plev(klon,klev+1),temp(klon,klev)
	35	real ustar(klon),snstable
	36	REAL zlev(klon,klev+1)
	37	REAL zlay(klon,klev)
	38	REAL u(klon,klev)
	39	REAL v(klon,klev)
	40	REAL teta(klon,klev)
	41	REAL cd(klon)
	42	REAL q2(klon,klev+1),q2s(klon,klev+1)
	43	REAL q2diag(klon,klev+1)
	44	REAL km(klon,klev+1)
	45	REAL kn(klon,klev+1)
	46	real sq(klon),sqz(klon),zz(klon,klev+1),zq,long0(klon)
	47
	48	integer l_mix,iii
	49	c.......................................................................
	50	c
	51	c nlay : nombre de couches
	52	c nlev : nombre de niveaux
	53	c ngrid : nombre de points de grille
	54	c unsdz : 1 sur l'epaisseur de couche
	55	c unsdzdec : 1 sur la distance entre le centre de la couche et le
	56	c centre de la couche inferieure
	57	c q : echelle de vitesse au bas de chaque couche
	58	c (valeur a la fin du pas de temps)
	59	c
	60	c.......................................................................
	61	INTEGER nlay,nlev,ngrid
	62	REAL unsdz(klon,klev)
	63	REAL unsdzdec(klon,klev+1)
	64	REAL q(klon,klev+1)
	65
	66	c.......................................................................
	67	c
	68	c kmpre : km au debut du pas de temps
	69	c qcstat : q : solution stationnaire du probleme couple
	70	c (valeur a la fin du pas de temps)
	71	c q2cstat : q2 : solution stationnaire du probleme couple
	72	c (valeur a la fin du pas de temps)
	73	c
	74	c.......................................................................
	75	REAL kmpre(klon,klev+1)
	76	REAL qcstat
	77	REAL q2cstat
	78	real sss,sssq
	79	c.......................................................................
	80	c
	81	c long : longueur de melange calculee selon Blackadar
	82	c
	83	c.......................................................................
	84	REAL long(klon,klev+1)
	85	c.......................................................................
	86	c
	87	c kmq3 : terme en q^3 dans le developpement de km
	88	c (valeur au debut du pas de temps)
	89	c kmcstat : valeur de km solution stationnaire du systeme {q2 ; du/dz}
	90	c (valeur a la fin du pas de temps)
	91	c knq3 : terme en q^3 dans le developpement de kn
	92	c mcstat : valeur de m solution stationnaire du systeme {q2 ; du/dz}
	93	c (valeur a la fin du pas de temps)
	94	c m2cstat : valeur de m2 solution stationnaire du systeme {q2 ; du/dz}
	95	c (valeur a la fin du pas de temps)
	96	c m : valeur a la fin du pas de temps
	97	c mpre : valeur au debut du pas de temps
	98	c m2 : valeur a la fin du pas de temps
	99	c n2 : valeur a la fin du pas de temps
	100	c
	101	c.......................................................................
	102	REAL kmq3
	103	REAL kmcstat
	104	REAL knq3
	105	REAL mcstat
	106	REAL m2cstat
	107	REAL m(klon,klev+1)
	108	REAL mpre(klon,klev+1)
	109	REAL m2(klon,klev+1)
	110	REAL n2(klon,klev+1)
	111	c.......................................................................
	112	c
	113	c gn : intermediaire pour les coefficients de stabilite
	114	c gnmin : borne inferieure de gn (-0.23 ou -0.28)
	115	c gnmax : borne superieure de gn (0.0233)
	116	c gninf : vrai si gn est en dessous de sa borne inferieure
	117	c gnsup : vrai si gn est en dessus de sa borne superieure
	118	c gm : drole d'objet bien utile
	119	c ri : nombre de Richardson
	120	c sn : coefficient de stabilite pour n
	121	c snq2 : premier terme du developement limite de sn en q2
	122	c sm : coefficient de stabilite pour m
	123	c smq2 : premier terme du developement limite de sm en q2
	124	c
	125	c.......................................................................
	126	REAL gn
	127	REAL gnmin
	128	REAL gnmax
	129	LOGICAL gninf
	130	LOGICAL gnsup
	131	REAL gm
	132	c REAL ri(klon,klev+1)
	133	REAL sn(klon,klev+1)
	134	REAL snq2(klon,klev+1)
	135	REAL sm(klon,klev+1)
	136	REAL smq2(klon,klev+1)
	137	c.......................................................................
	138	c
	139	c kappa : consatnte de Von Karman (0.4)
	140	c long00 : longueur de reference pour le calcul de long (160)
	141	c a1,a2,b1,b2,c1 : constantes d'origine pour les coefficients
	142	c de stabilite (0.92/0.74/16.6/10.1/0.08)
	143	c cn1,cn2 : constantes pour sn
	144	c cm1,cm2,cm3,cm4 : constantes pour sm
	145	c
	146	c.......................................................................
	147	REAL kappa
	148	REAL long00
	149	REAL a1,a2,b1,b2,c1
	150	REAL cn1,cn2
	151	REAL cm1,cm2,cm3,cm4
	152	c.......................................................................
	153	c
	154	c termq : termes en $q$ dans l'equation de q2
	155	c termq3 : termes en $q^3$ dans l'equation de q2
	156	c termqm2 : termes en $q*m^2$ dans l'equation de q2
	157	c termq3m2 : termes en $q^3*m^2$ dans l'equation de q2
	158	c
	159	c.......................................................................
	160	REAL termq
	161	REAL termq3
	162	REAL termqm2
	163	REAL termq3m2
	164	c.......................................................................
	165	c
	166	c q2min : borne inferieure de q2
	167	c q2max : borne superieure de q2
	168	c
	169	c.......................................................................
	170	REAL q2min
	171	REAL q2max
	172	c.......................................................................
	173	c knmin : borne inferieure de kn
	174	c kmmin : borne inferieure de km
	175	c.......................................................................
	176	REAL knmin
	177	REAL kmmin
	178	c.......................................................................
	179	INTEGER ilay,ilev,igrid
	180	REAL tmp1,tmp2
	181	c.......................................................................
	182	PARAMETER (kappa=0.4E+0)
	183	PARAMETER (long00=160.E+0)
	184	c PARAMETER (gnmin=-10.E+0)
	185	PARAMETER (gnmin=-0.28)
	186	PARAMETER (gnmax=0.0233E+0)
	187	PARAMETER (a1=0.92E+0)
	188	PARAMETER (a2=0.74E+0)
	189	PARAMETER (b1=16.6E+0)
	190	PARAMETER (b2=10.1E+0)
	191	PARAMETER (c1=0.08E+0)
	192	PARAMETER (knmin=1.E-5)
	193	PARAMETER (kmmin=1.E-5)
	194	PARAMETER (q2min=1.e-5)
	195	PARAMETER (q2max=1.E+2)
	196	c
	197	PARAMETER (
	198	& cn1=a2(1.E+0 -6.E+0 a1/b1)
	199	& )
	200	PARAMETER (
	201	& cn2=-3.E+0 a2(6.E+0 *a1+b2)
	202	& )
	203	PARAMETER (
	204	& cm1=a1(1.E+0 -3.E+0 c1-6.E+0 *a1/b1)
	205	& )
	206	PARAMETER (
	207	& cm2=a1(-3.E+0 a2((b2-3.E+0 a2)(1.E+0 -6.E+0 a1/b1)
	208	& -3.E+0 c1(b2+6.E+0 *a1)))
	209	& )
	210	PARAMETER (
	211	& cm3=-3.E+0 a2(6.E+0 *a1+b2)
	212	& )
	213	PARAMETER (
	214	& cm4=-9.E+0 a1a2
	215	& )
	216
	217	logical first
	218	save first
	219	data first/.true./
[101]	220
	221	nlay=klev
	222	nlev=klev+1
[3]	223	c.......................................................................
	224	c traitment des valeur de q2 en entree
	225	c.......................................................................
	226	c
	227	c Initialisation de q2
	228
	229	call yamada(ngrid,dt,g,rconst,plev,temp
	230	s ,zlev,zlay,u,v,teta,cd,q2diag,km,kn,ustar
	231	s ,l_mix)
	232	if (first.and.1.eq.1) then
	233	first=.false.
	234	q2=q2diag
	235	endif
	236
	237	DO ilev=1,nlev
	238	DO igrid=1,ngrid
	239	q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
	240	q(igrid,ilev)=sqrt(q2(igrid,ilev))
	241	ENDDO
	242	ENDDO
	243	c
	244	DO igrid=1,ngrid
	245	tmp1=cd(igrid)(u(igrid,1)2+v(igrid,1)*2)
	246	q2(igrid,1)=b1*(2.E+0/3.E+0)tmp1
	247	q2(igrid,1)=amax1(q2(igrid,1),q2min)
	248	q(igrid,1)=sqrt(q2(igrid,1))
	249	ENDDO
	250	c
	251	c.......................................................................
	252	c les increments verticaux
	253	c.......................................................................
	254	c
	255	c!!!!! allerte !!!!!c
	256	c!!!!! zlev n'est pas declare a nlev !!!!!c
	257	c!!!!! ---->
	258	DO igrid=1,ngrid
	259	zlev(igrid,nlev)=zlay(igrid,nlay)
	260	& +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
	261	ENDDO
	262	c!!!!! <----
	263	c!!!!! allerte !!!!!c
	264	c
	265	DO ilay=1,nlay
	266	DO igrid=1,ngrid
	267	unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
	268	ENDDO
	269	ENDDO
	270	DO igrid=1,ngrid
	271	unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
	272	ENDDO
	273	DO ilay=2,nlay
	274	DO igrid=1,ngrid
	275	unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
	276	ENDDO
	277	ENDDO
	278	DO igrid=1,ngrid
	279	unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
	280	ENDDO
	281	c
	282	c.......................................................................
	283	c le cisaillement et le gradient de temperature
	284	c.......................................................................
	285	c
	286	DO igrid=1,ngrid
	287	m2(igrid,1)=(unsdzdec(igrid,1)
	288	& u(igrid,1))*2
	289	& +(unsdzdec(igrid,1)
	290	& v(igrid,1))*2
	291	m(igrid,1)=sqrt(m2(igrid,1))
	292	mpre(igrid,1)=m(igrid,1)
	293	ENDDO
	294	c
	295	c-----------------------------------------------------------------------
	296	DO ilev=2,nlev-1
	297	DO igrid=1,ngrid
	298	c-----------------------------------------------------------------------
	299	c
	300	n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
	301	& *(teta(igrid,ilev)-teta(igrid,ilev-1))
	302	& /(teta(igrid,ilev)+teta(igrid,ilev-1)) *2.E+0
	303	c n2(igrid,ilev)=0.
	304	c
	305	c --->
	306	c on ne sais traiter que les cas stratifies. et l'ajustement
	307	c convectif est cense faire en sorte que seul des configurations
	308	c stratifiees soient rencontrees en entree de cette routine.
	309	c mais, bon ... on sait jamais (meme on sait que n2 prends
	310	c quelques valeurs negatives ... parfois) alors :
	311	c<---
	312	c
	313	IF (n2(igrid,ilev).lt.0.E+0) THEN
	314	n2(igrid,ilev)=0.E+0
	315	ENDIF
	316	c
	317	m2(igrid,ilev)=(unsdzdec(igrid,ilev)
	318	& (u(igrid,ilev)-u(igrid,ilev-1)))*2
	319	& +(unsdzdec(igrid,ilev)
	320	& (v(igrid,ilev)-v(igrid,ilev-1)))*2
	321	m(igrid,ilev)=sqrt(m2(igrid,ilev))
	322	mpre(igrid,ilev)=m(igrid,ilev)
	323	c
	324	c-----------------------------------------------------------------------
	325	ENDDO
	326	ENDDO
	327	c-----------------------------------------------------------------------
	328	c
	329	DO igrid=1,ngrid
	330	m2(igrid,nlev)=m2(igrid,nlev-1)
	331	m(igrid,nlev)=m(igrid,nlev-1)
	332	mpre(igrid,nlev)=m(igrid,nlev)
	333	ENDDO
	334	c
	335	c.......................................................................
	336	c calcul des fonctions de stabilite
	337	c.......................................................................
	338	c
	339	if (l_mix.eq.4) then
	340	DO igrid=1,ngrid
	341	sqz(igrid)=1.e-10
	342	sq(igrid)=1.e-10
	343	ENDDO
	344	do ilev=2,nlev-1
	345	DO igrid=1,ngrid
	346	zq=sqrt(q2(igrid,ilev))
	347	sqz(igrid)
	348	. =sqz(igrid)+zq*zlev(igrid,ilev)
	349	. *(zlay(igrid,ilev)-zlay(igrid,ilev-1))
	350	sq(igrid)=sq(igrid)+zq*(zlay(igrid,ilev)-zlay(igrid,ilev-1))
	351	ENDDO
	352	enddo
	353	DO igrid=1,ngrid
	354	long0(igrid)=0.2*sqz(igrid)/sq(igrid)
	355	ENDDO
	356	else if (l_mix.eq.3) then
	357	long0(igrid)=long00
	358	endif
	359
	360	c (abd 5 2) print*,'LONG0=',long0
	361
	362	c-----------------------------------------------------------------------
	363	DO ilev=2,nlev-1
	364	DO igrid=1,ngrid
	365	c-----------------------------------------------------------------------
	366	c
	367	tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
	368	if (l_mix.ge.10) then
	369	long(igrid,ilev)=l_mix
	370	else
	371	long(igrid,ilev)=tmp1/(1.E+0 + tmp1/long0(igrid))
	372	endif
	373	long(igrid,ilev)=max(min(long(igrid,ilev)
	374	s ,0.5*sqrt(q2(igrid,ilev))/sqrt(max(n2(igrid,ilev),1.e-10)))
	375	s ,5.)
	376
	377	gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
	378	& * n2(igrid,ilev)
	379	gm=long(igrid,ilev)**2 / q2(igrid,ilev)
	380	& * m2(igrid,ilev)
	381	c
	382	gninf=.false.
	383	gnsup=.false.
	384	long(igrid,ilev)=long(igrid,ilev)
	385	long(igrid,ilev)=long(igrid,ilev)
	386	c
	387	IF (gn.lt.gnmin) THEN
	388	gninf=.true.
	389	gn=gnmin
	390	ENDIF
	391	c
	392	IF (gn.gt.gnmax) THEN
	393	gnsup=.true.
	394	gn=gnmax
	395	ENDIF
	396	c
	397	sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
	398	sm(igrid,ilev)=
	399	& (cm1+cm2*gn)
	400	& /( (1.E+0 +cm3*gn)
	401	& (1.E+0 +cm4gn) )
	402	c
	403	IF ((gninf).or.(gnsup)) THEN
	404	snq2(igrid,ilev)=0.E+0
	405	smq2(igrid,ilev)=0.E+0
	406	ELSE
	407	snq2(igrid,ilev)=
	408	& -gn
	409	& (-cn1cn2/(1.E+0 +cn2gn)*2 )
	410	smq2(igrid,ilev)=
	411	& -gn
	412	& ( cm2(1.E+0 +cm3*gn)
	413	& (1.E+0 +cm4gn)
	414	& -( cm3(1.E+0 +cm4gn)
	415	& +cm4(1.E+0 +cm3gn) )
	416	& (cm1+cm2gn) )
	417	& /( (1.E+0 +cm3*gn)
	418	& (1.E+0 +cm4gn) )**2
	419	ENDIF
	420	c
	421	c abd
	422	c if(ilev.le.57.and.ilev.ge.37) then
	423	c print*,'L=',ilev,' GN=',gn,' SM=',sm(igrid,ilev)
	424	c endif
	425	c --->
	426	c la decomposition de Taylor en q2 n'a de sens que
	427	c dans les cas stratifies ou sn et sm sont quasi
	428	c proportionnels a q2. ailleurs on laisse le meme
	429	c algorithme car l'ajustement convectif fait le travail.
	430	c mais c'est delirant quand sn et snq2 n'ont pas le meme
	431	c signe : dans ces cas, on ne fait pas la decomposition.
	432	c<---
	433	c
	434	IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.E+0)
	435	& snq2(igrid,ilev)=0.E+0
	436	IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.E+0)
	437	& smq2(igrid,ilev)=0.E+0
	438	c
	439	C Correction pour les couches stables.
	440	C Schema repris de JHoltzlag Boville, lui meme venant de...
	441
	442	if (1.eq.1) then
	443	snstable=1.-zlev(igrid,ilev)
	444	s /(700.*max(ustar(igrid),0.0001))
	445	snstable=1.-zlev(igrid,ilev)/400.
	446	snstable=max(snstable,0.)
	447	snstable=snstable*snstable
	448
	449	c abde print*,'SN ',ilev,sn(1,ilev),snstable
	450	if (sn(igrid,ilev).lt.snstable) then
	451	sn(igrid,ilev)=snstable
	452	snq2(igrid,ilev)=0.
	453	endif
	454
	455	if (sm(igrid,ilev).lt.snstable) then
	456	sm(igrid,ilev)=snstable
	457	smq2(igrid,ilev)=0.
	458	endif
	459
	460	endif
	461
	462	c sn : coefficient de stabilite pour n
	463	c snq2 : premier terme du developement limite de sn en q2
	464	c-----------------------------------------------------------------------
	465	ENDDO
	466	ENDDO
	467	c-----------------------------------------------------------------------
	468	c
	469	c.......................................................................
	470	c calcul de km et kn au debut du pas de temps
	471	c.......................................................................
	472	c
	473	DO igrid=1,ngrid
	474	kn(igrid,1)=knmin
	475	km(igrid,1)=kmmin
	476	kmpre(igrid,1)=km(igrid,1)
	477	ENDDO
	478	c
	479	c-----------------------------------------------------------------------
	480	DO ilev=2,nlev-1
	481	DO igrid=1,ngrid
	482	c-----------------------------------------------------------------------
	483	c
	484	kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	485	& *sn(igrid,ilev)
	486	km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	487	& *sm(igrid,ilev)
	488	kmpre(igrid,ilev)=km(igrid,ilev)
	489	c
	490	c-----------------------------------------------------------------------
	491	ENDDO
	492	ENDDO
	493	c-----------------------------------------------------------------------
	494	c
	495	DO igrid=1,ngrid
	496	kn(igrid,nlev)=kn(igrid,nlev-1)
	497	km(igrid,nlev)=km(igrid,nlev-1)
	498	kmpre(igrid,nlev)=km(igrid,nlev)
	499	ENDDO
	500	c
	501	c.......................................................................
	502	c boucle sur les niveaux 2 a nlev-1
	503	c.......................................................................
	504	c
	505	c---->
	506	DO 10001 ilev=2,nlev-1
	507	c---->
	508	DO 10002 igrid=1,ngrid
	509	c
	510	c.......................................................................
	511	c
	512	c calcul des termes sources et puits de l'equation de q2
	513	c ------------------------------------------------------
	514	c
	515	knq3=kn(igrid,ilev)*snq2(igrid,ilev)
	516	& /sn(igrid,ilev)
	517	kmq3=km(igrid,ilev)*smq2(igrid,ilev)
	518	& /sm(igrid,ilev)
	519	c
	520	termq=0.E+0
	521	termq3=0.E+0
	522	termqm2=0.E+0
	523	termq3m2=0.E+0
	524	c
	525	tmp1=dt2.E+0 km(igrid,ilev)*m2(igrid,ilev)
	526	tmp2=dt2.E+0 kmq3*m2(igrid,ilev)
	527	termqm2=termqm2
	528	& +dt2.E+0 km(igrid,ilev)*m2(igrid,ilev)
	529	& -dt2.E+0 kmq3*m2(igrid,ilev)
	530	termq3m2=termq3m2
	531	& +dt2.E+0 kmq3*m2(igrid,ilev)
	532	c
	533	termq=termq
	534	& -dt2.E+0 kn(igrid,ilev)*n2(igrid,ilev)
	535	& +dt2.E+0 knq3*n2(igrid,ilev)
	536	termq3=termq3
	537	& -dt2.E+0 knq3*n2(igrid,ilev)
	538	c
	539	termq3=termq3
	540	& -dt2.E+0 q(igrid,ilev)*3 / (b1long(igrid,ilev))
	541	c
	542	c.......................................................................
	543	c
	544	c resolution stationnaire couplee avec le gradient de vitesse local
	545	c -----------------------------------------------------------------
	546	c
	547	c -----{on cherche le cisaillement qui annule l'equation de q^2
	548	c supposee en q3}
	549	c
	550	tmp1=termq+termq3
	551	tmp2=termqm2+termq3m2
	552	m2cstat=m2(igrid,ilev)
	553	& -(tmp1+tmp2)/(dt2.E+0km(igrid,ilev))
	554	mcstat=sqrt(m2cstat)
	555
	556	c abde print*,'M2 L=',ilev,mpre(igrid,ilev),mcstat
	557	c
	558	c -----{puis on ecrit la valeur de q qui annule l'equation de m
	559	c supposee en q3}
	560	c
	561	IF (ilev.eq.2) THEN
	562	kmcstat=1.E+0 / mcstat
	563	& ( unsdz(igrid,ilev)kmpre(igrid,ilev+1)
	564	& *mpre(igrid,ilev+1)
	565	& +unsdz(igrid,ilev-1)
	566	& *cd(igrid)
	567	& ( sqrt(u(igrid,3)2+v(igrid,3)*2)
	568	& -mcstat/unsdzdec(igrid,ilev)
	569	& -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
	570	& /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
	571	ELSE
	572	kmcstat=1.E+0 / mcstat
	573	& ( unsdz(igrid,ilev)kmpre(igrid,ilev+1)
	574	& *mpre(igrid,ilev+1)
	575	& +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
	576	& *mpre(igrid,ilev-1) )
	577	& /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
	578	ENDIF
	579	tmp2=kmcstat
	580	& /( sm(igrid,ilev)/q2(igrid,ilev) )
	581	& /long(igrid,ilev)
	582	qcstat=tmp2**(1.E+0/3.E+0)
	583	q2cstat=qcstat**2
	584	c
	585	c.......................................................................
	586	c
	587	c choix de la solution finale
	588	c ---------------------------
	589	c
	590	q(igrid,ilev)=qcstat
	591	q2(igrid,ilev)=q2cstat
	592	m(igrid,ilev)=mcstat
	593	c abd if(ilev.le.57.and.ilev.ge.37) then
	594	c print*,'L=',ilev,' M2=',m2(igrid,ilev),m2cstat,
	595	c s 'N2=',n2(igrid,ilev)
	596	c abd endif
	597	m2(igrid,ilev)=m2cstat
	598	c
	599	c --->
	600	c pour des raisons simples q2 est minore
	601	c<---
	602	c
	603	IF (q2(igrid,ilev).lt.q2min) THEN
	604	q2(igrid,ilev)=q2min
	605	q(igrid,ilev)=sqrt(q2min)
	606	ENDIF
	607	c
	608	c.......................................................................
	609	c
	610	c calcul final de kn et km
	611	c ------------------------
	612	c
	613	gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
	614	& * n2(igrid,ilev)
	615	IF (gn.lt.gnmin) gn=gnmin
	616	IF (gn.gt.gnmax) gn=gnmax
	617	sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
	618	sm(igrid,ilev)=
	619	& (cm1+cm2*gn)
	620	& /( (1.E+0 +cm3gn)(1.E+0 +cm4*gn) )
	621	kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	622	& *sn(igrid,ilev)
	623	km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	624	& *sm(igrid,ilev)
	625	c abd
	626	c if(ilev.le.57.and.ilev.ge.37) then
	627	c print*,'L=',ilev,' GN=',gn,' SM=',sm(igrid,ilev)
	628	c endif
	629	c
	630	c.......................................................................
	631	c
	632	10002 CONTINUE
	633	c
	634	10001 CONTINUE
	635	c
	636	c.......................................................................
	637	c
	638	c
	639	DO igrid=1,ngrid
	640	kn(igrid,1)=knmin
	641	km(igrid,1)=kmmin
	642	c kn(igrid,1)=cd(igrid)
	643	c km(igrid,1)=cd(igrid)
	644	q2(igrid,nlev)=q2(igrid,nlev-1)
	645	q(igrid,nlev)=q(igrid,nlev-1)
	646	kn(igrid,nlev)=kn(igrid,nlev-1)
	647	km(igrid,nlev)=km(igrid,nlev-1)
	648	ENDDO
	649	c
	650	c CALCUL DE LA DIFFUSION VERTICALE DE Q2
	651	if (1.eq.1) then
	652
	653	do ilev=2,klev-1
	654	sss=sss+plev(1,ilev-1)-plev(1,ilev+1)
	655	sssq=sssq+(plev(1,ilev-1)-plev(1,ilev+1))*q2(1,ilev)
	656	enddo
	657	c print*,'Q2moy avant',sssq/sss
	658	c print*,'Q2q20 ',(q2(1,ilev),ilev=1,10)
	659	c print*,'Q2km0 ',(km(1,ilev),ilev=1,10)
	660	c ! C'est quoi ca qu'etait dans l'original???
	661	c do igrid=1,ngrid
	662	c q2(igrid,1)=10.
	663	c enddo
	664	c q2s=q2
	665	c do iii=1,10
	666	c call vdif_q2(dt,g,rconst,plev,temp,km,q2)
	667	c do ilev=1,klev+1
	668	c write(iii+49,*) q2(1,ilev),zlev(1,ilev)
	669	c enddo
	670	c enddo
	671	c stop
	672	c do ilev=1,klev
	673	c print*,zlev(1,ilev),q2s(1,ilev),q2(1,ilev)
	674	c enddo
	675	c q2s=q2-q2s
	676	c do ilev=1,klev
	677	c print*,q2s(1,ilev),zlev(1,ilev)
	678	c enddo
	679	do ilev=2,klev-1
	680	sss=sss+plev(1,ilev-1)-plev(1,ilev+1)
	681	sssq=sssq+(plev(1,ilev-1)-plev(1,ilev+1))*q2(1,ilev)
	682	enddo
	683	print*,'Q2moy apres',sssq/sss
	684	c
	685	c
	686	do ilev=1,nlev
	687	do igrid=1,ngrid
	688	q2(igrid,ilev)=max(q2(igrid,ilev),q2min)
	689	q(igrid,ilev)=sqrt(q2(igrid,ilev))
	690
	691	c.......................................................................
	692	c
	693	c calcul final de kn et km
	694	c ------------------------
	695	c
	696	gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
	697	& * n2(igrid,ilev)
	698	IF (gn.lt.gnmin) gn=gnmin
	699	IF (gn.gt.gnmax) gn=gnmax
	700	sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
	701	sm(igrid,ilev)=
	702	& (cm1+cm2*gn)
	703	& /( (1.E+0 +cm3gn)(1.E+0 +cm4*gn) )
	704	C Correction pour les couches stables.
	705	C Schema repris de JHoltzlag Boville, lui meme venant de...
	706
	707	if (1.eq.1) then
	708	snstable=1.-zlev(igrid,ilev)
	709	s /(700.*max(ustar(igrid),0.0001))
	710	snstable=1.-zlev(igrid,ilev)/400.
	711	snstable=max(snstable,0.)
	712	snstable=snstable*snstable
	713
	714	c abde print*,'SN ',ilev,sn(1,ilev),snstable
	715	if (sn(igrid,ilev).lt.snstable) then
	716	sn(igrid,ilev)=snstable
	717	snq2(igrid,ilev)=0.
	718	endif
	719
	720	if (sm(igrid,ilev).lt.snstable) then
	721	sm(igrid,ilev)=snstable
	722	smq2(igrid,ilev)=0.
	723	endif
	724
	725	endif
	726
	727	c sn : coefficient de stabilite pour n
	728	kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	729	& *sn(igrid,ilev)
	730	km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
	731	c
	732	enddo
	733	enddo
	734	c print*,'Q2km1 ',(km(1,ilev),ilev=1,10)
	735
	736	endif
	737
	738	RETURN
	739	END

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

Download in other formats:

Original Format