Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

vdif_kcay.F90 @ 1531

Last change on this file since 1531 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: 29.5 KB

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

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

Download in other formats:

Original Format