Context Navigation

vdif_kcay.F @ 1198

Last change on this file since 1198 was 101, checked in by slebonnois, 14 years ago
SL: modifications pour arriver a compiler le gcm VENUS ! Ca marche ! A noter: modifs de makelmdz
File size: 25.1 KB

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

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

Download in other formats:

Original Format