Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

vdif_kc.F @ 38

Last change on this file since 38 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 19.8 KB

Line
1	SUBROUTINE vdif_kc(dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
2	IMPLICIT NONE
3	c.......................................................................
4	#include "dimensions.h"
5	#include "dimphys.h"
6	c.......................................................................
7	c
8	c dt : pas de temps
9	c g : g
10	c zlev : altitude a chaque niveau (interface inferieure de la couche
11	c de meme indice)
12	c zlay : altitude au centre de chaque couche
13	c u,v : vitesse au centre de chaque couche
14	c (en entree : la valeur au debut du pas de temps)
15	c teta : temperature potentielle au centre de chaque couche
16	c (en entree : la valeur au debut du pas de temps)
17	c cd : cdrag
18	c (en entree : la valeur au debut du pas de temps)
19	c q2 : $q^2$ au bas de chaque couche
20	c (en entree : la valeur au debut du pas de temps)
21	c (en sortie : la valeur a la fin du pas de temps)
22	c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
23	c couche)
24	c (en sortie : la valeur a la fin du pas de temps)
25	c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
26	c (en sortie : la valeur a la fin du pas de temps)
27	c
28	c.......................................................................
29	REAL dt,g
30	REAL zlev(ngridmx,nlayermx+1)
31	REAL zlay(ngridmx,nlayermx)
32	REAL u(ngridmx,nlayermx)
33	REAL v(ngridmx,nlayermx)
34	REAL teta(ngridmx,nlayermx)
35	REAL cd(ngridmx)
36	REAL q2(ngridmx,nlayermx+1)
37	REAL km(ngridmx,nlayermx+1)
38	REAL kn(ngridmx,nlayermx+1)
39	c.......................................................................
40	c
41	c nlay : nombre de couches
42	c nlev : nombre de niveaux
43	c ngrid : nombre de points de grille
44	c unsdz : 1 sur l'epaisseur de couche
45	c unsdzdec : 1 sur la distance entre le centre de la couche et le
46	c centre de la couche inferieure
47	c q : echelle de vitesse au bas de chaque couche
48	c (valeur a la fin du pas de temps)
49	c
50	c.......................................................................
51	INTEGER nlay,nlev,ngrid
52	REAL unsdz(ngridmx,nlayermx)
53	REAL unsdzdec(ngridmx,nlayermx+1)
54	REAL q(ngridmx,nlayermx+1)
55	c.......................................................................
56	c
57	c kmpre : km au debut du pas de temps
58	c qcstat : q : solution stationnaire du probleme couple
59	c (valeur a la fin du pas de temps)
60	c q2cstat : q2 : solution stationnaire du probleme couple
61	c (valeur a la fin du pas de temps)
62	c
63	c.......................................................................
64	REAL kmpre(ngridmx,nlayermx+1)
65	REAL qcstat
66	REAL q2cstat
67	c.......................................................................
68	c
69	c long : longueur de melange calculee selon Blackadar
70	c
71	c.......................................................................
72	REAL long(ngridmx,nlayermx+1)
73	c.......................................................................
74	c
75	c kmq3 : terme en q^3 dans le developpement de km
76	c (valeur au debut du pas de temps)
77	c kmcstat : valeur de km solution stationnaire du systeme {q2 ; du/dz}
78	c (valeur a la fin du pas de temps)
79	c knq3 : terme en q^3 dans le developpement de kn
80	c mcstat : valeur de m solution stationnaire du systeme {q2 ; du/dz}
81	c (valeur a la fin du pas de temps)
82	c m2cstat : valeur de m2 solution stationnaire du systeme {q2 ; du/dz}
83	c (valeur a la fin du pas de temps)
84	c m : valeur a la fin du pas de temps
85	c mpre : valeur au debut du pas de temps
86	c m2 : valeur a la fin du pas de temps
87	c n2 : valeur a la fin du pas de temps
88	c
89	c.......................................................................
90	REAL kmq3
91	REAL kmcstat
92	REAL knq3
93	REAL mcstat
94	REAL m2cstat
95	REAL m(ngridmx,nlayermx+1)
96	REAL mpre(ngridmx,nlayermx+1)
97	REAL m2(ngridmx,nlayermx+1)
98	REAL n2(ngridmx,nlayermx+1)
99	c.......................................................................
100	c
101	c gn : intermediaire pour les coefficients de stabilite
102	c gnmin : borne inferieure de gn (-0.23 ou -0.28)
103	c gnmax : borne superieure de gn (0.0233)
104	c gninf : vrai si gn est en dessous de sa borne inferieure
105	c gnsup : vrai si gn est en dessus de sa borne superieure
106	c gm : drole d'objet bien utile
107	c ri : nombre de Richardson
108	c sn : coefficient de stabilite pour n
109	c snq2 : premier terme du developement limite de sn en q2
110	c sm : coefficient de stabilite pour m
111	c smq2 : premier terme du developement limite de sm en q2
112	c
113	c.......................................................................
114	REAL gn
115	REAL gnmin
116	REAL gnmax
117	LOGICAL gninf
118	LOGICAL gnsup
119	REAL gm
120	c REAL ri(ngridmx,nlayermx+1)
121	REAL sn(ngridmx,nlayermx+1)
122	REAL snq2(ngridmx,nlayermx+1)
123	REAL sm(ngridmx,nlayermx+1)
124	REAL smq2(ngridmx,nlayermx+1)
125	c.......................................................................
126	c
127	c kappa : consatnte de Von Karman (0.4)
128	c long0 : longueur de reference pour le calcul de long (160)
129	c a1,a2,b1,b2,c1 : constantes d'origine pour les coefficients
130	c de stabilite (0.92/0.74/16.6/10.1/0.08)
131	c cn1,cn2 : constantes pour sn
132	c cm1,cm2,cm3,cm4 : constantes pour sm
133	c
134	c.......................................................................
135	REAL kappa
136	REAL long0
137	REAL a1,a2,b1,b2,c1
138	REAL cn1,cn2
139	REAL cm1,cm2,cm3,cm4
140	c.......................................................................
141	c
142	c termq : termes en $q$ dans l'equation de q2
143	c termq3 : termes en $q^3$ dans l'equation de q2
144	c termqm2 : termes en $q*m^2$ dans l'equation de q2
145	c termq3m2 : termes en $q^3*m^2$ dans l'equation de q2
146	c
147	c.......................................................................
148	REAL termq
149	REAL termq3
150	REAL termqm2
151	REAL termq3m2
152	c.......................................................................
153	c
154	c q2min : borne inferieure de q2
155	c q2max : borne superieure de q2
156	c
157	c.......................................................................
158	REAL q2min
159	REAL q2max
160	c.......................................................................
161	c knmin : borne inferieure de kn
162	c kmmin : borne inferieure de km
163	c.......................................................................
164	REAL knmin
165	REAL kmmin
166	c.......................................................................
167	INTEGER ilay,ilev,igrid
168	REAL tmp1,tmp2
169	c.......................................................................
170	PARAMETER (kappa=0.4E+0)
171	PARAMETER (long0=160.E+0)
172	PARAMETER (gnmin=-10.E+0)
173	PARAMETER (gnmax=0.0233E+0)
174	PARAMETER (a1=0.92E+0)
175	PARAMETER (a2=0.74E+0)
176	PARAMETER (b1=16.6E+0)
177	PARAMETER (b2=10.1E+0)
178	PARAMETER (c1=0.08E+0)
179	PARAMETER (knmin=1.E-5)
180	PARAMETER (kmmin=1.E-5)
181	PARAMETER (q2min=1.E-3)
182	PARAMETER (q2max=1.E+2)
183	PARAMETER (nlay=nlayermx)
184	PARAMETER (nlev=nlayermx+1)
185	PARAMETER (ngrid=ngridmx)
186	c
187	PARAMETER (
188	& cn1=a2(1.E+0 -6.E+0 a1/b1)
189	& )
190	PARAMETER (
191	& cn2=-3.E+0 a2(6.E+0 *a1+b2)
192	& )
193	PARAMETER (
194	& cm1=a1(1.E+0 -3.E+0 c1-6.E+0 *a1/b1)
195	& )
196	PARAMETER (
197	& cm2=a1(-3.E+0 a2((b2-3.E+0 a2)(1.E+0 -6.E+0 a1/b1)
198	& -3.E+0 c1(b2+6.E+0 *a1)))
199	& )
200	PARAMETER (
201	& cm3=-3.E+0 a2(6.E+0 *a1+b2)
202	& )
203	PARAMETER (
204	& cm4=-9.E+0 a1a2
205	& )
206	c.......................................................................
207	c traitment des valeur de q2 en entree
208	c.......................................................................
209	c
210	DO ilev=1,nlev
211	DO igrid=1,ngrid
212	q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
213	q(igrid,ilev)=sqrt(q2(igrid,ilev))
214	ENDDO
215	ENDDO
216	c
217	DO igrid=1,ngrid
218	tmp1=cd(igrid)(u(igrid,1)2+v(igrid,1)*2)
219	q2(igrid,1)=b1*(2.E+0/3.E+0)tmp1
220	q2(igrid,1)=amax1(q2(igrid,1),q2min)
221	q(igrid,1)=sqrt(q2(igrid,1))
222	ENDDO
223	c
224	c.......................................................................
225	c les increments verticaux
226	c.......................................................................
227	c
228	c!!!!! allerte !!!!!c
229	c!!!!! zlev n'est pas declare a nlev !!!!!c
230	c!!!!! ---->
231	c DO igrid=1,ngrid
232	c zlev(igrid,nlev)=zlay(igrid,nlay)
233	c & +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
234	c ENDDO
235	c!!!!! <----
236	c!!!!! allerte !!!!!c
237	c
238	DO ilay=1,nlay
239	DO igrid=1,ngrid
240	unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
241	ENDDO
242	ENDDO
243	DO igrid=1,ngrid
244	unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
245	ENDDO
246	DO ilay=2,nlay
247	DO igrid=1,ngrid
248	unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
249	ENDDO
250	ENDDO
251	DO igrid=1,ngrid
252	unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
253	ENDDO
254	c
255	c.......................................................................
256	c le cisaillement et le gradient de temperature
257	c.......................................................................
258	c
259	DO igrid=1,ngrid
260	m2(igrid,1)=(unsdzdec(igrid,1)
261	& u(igrid,1))*2
262	& +(unsdzdec(igrid,1)
263	& v(igrid,1))*2
264	m(igrid,1)=sqrt(m2(igrid,1))
265	mpre(igrid,1)=m(igrid,1)
266	ENDDO
267	c
268	c-----------------------------------------------------------------------
269	DO ilev=2,nlev-1
270	DO igrid=1,ngrid
271	c-----------------------------------------------------------------------
272	c
273	n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
274	& *(teta(igrid,ilev)-teta(igrid,ilev-1))
275	& /(teta(igrid,ilev)+teta(igrid,ilev-1)) *2.E+0
276	c
277	c --->
278	c on ne sais traiter que les cas stratifies. et l'ajustement
279	c convectif est cense faire en sorte que seul des configurations
280	c stratifiees soient rencontrees en entree de cette routine.
281	c mais, bon ... on sait jamais (meme on sait que n2 prends
282	c quelques valeurs negatives ... parfois) alors :
283	c<---
284	c
285	IF (n2(igrid,ilev).lt.0.E+0) THEN
286	n2(igrid,ilev)=0.E+0
287	ENDIF
288	c
289	m2(igrid,ilev)=(unsdzdec(igrid,ilev)
290	& (u(igrid,ilev)-u(igrid,ilev-1)))*2
291	& +(unsdzdec(igrid,ilev)
292	& (v(igrid,ilev)-v(igrid,ilev-1)))*2
293	m(igrid,ilev)=sqrt(m2(igrid,ilev))
294	mpre(igrid,ilev)=m(igrid,ilev)
295	c
296	c-----------------------------------------------------------------------
297	ENDDO
298	ENDDO
299	c-----------------------------------------------------------------------
300	c
301	DO igrid=1,ngrid
302	m2(igrid,nlev)=m2(igrid,nlev-1)
303	m(igrid,nlev)=m(igrid,nlev-1)
304	mpre(igrid,nlev)=m(igrid,nlev)
305	ENDDO
306	c
307	c.......................................................................
308	c calcul des fonctions de stabilite
309	c.......................................................................
310	c
311	c-----------------------------------------------------------------------
312	DO ilev=2,nlev-1
313	DO igrid=1,ngrid
314	c-----------------------------------------------------------------------
315	c
316	tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
317	long(igrid,ilev)=tmp1/(1.E+0 + tmp1/long0)
318	gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
319	& * n2(igrid,ilev)
320	gm=long(igrid,ilev)**2 / q2(igrid,ilev)
321	& * m2(igrid,ilev)
322	c
323	gninf=.false.
324	gnsup=.false.
325	long(igrid,ilev)=long(igrid,ilev)
326	long(igrid,ilev)=long(igrid,ilev)
327	c
328	IF (gn.lt.gnmin) THEN
329	gninf=.true.
330	gn=gnmin
331	ENDIF
332	c
333	IF (gn.gt.gnmax) THEN
334	gnsup=.true.
335	gn=gnmax
336	ENDIF
337	c
338	sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
339	sm(igrid,ilev)=
340	& (cm1+cm2*gn)
341	& /( (1.E+0 +cm3*gn)
342	& (1.E+0 +cm4gn) )
343	c
344	IF ((gninf).or.(gnsup)) THEN
345	snq2(igrid,ilev)=0.E+0
346	smq2(igrid,ilev)=0.E+0
347	ELSE
348	snq2(igrid,ilev)=
349	& -gn
350	& (-cn1cn2/(1.E+0 +cn2gn)*2 )
351	smq2(igrid,ilev)=
352	& -gn
353	& ( cm2(1.E+0 +cm3*gn)
354	& (1.E+0 +cm4gn)
355	& -( cm3(1.E+0 +cm4gn)
356	& +cm4(1.E+0 +cm3gn) )
357	& (cm1+cm2gn) )
358	& /( (1.E+0 +cm3*gn)
359	& (1.E+0 +cm4gn) )**2
360	ENDIF
361	c
362	c --->
363	c la decomposition de Taylor en q2 n'a de sens que
364	c dans les cas stratifies ou sn et sm sont quasi
365	c proportionnels a q2. ailleurs on laisse le meme
366	c algorithme car l'ajustement convectif fait le travail.
367	c mais c'est delirant quand sn et snq2 n'ont pas le meme
368	c signe : dans ces cas, on ne fait pas la decomposition.
369	c<---
370	c
371	IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.E+0)
372	& snq2(igrid,ilev)=0.E+0
373	IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.E+0)
374	& smq2(igrid,ilev)=0.E+0
375	c
376	c-----------------------------------------------------------------------
377	ENDDO
378	ENDDO
379	c-----------------------------------------------------------------------
380	c
381	c.......................................................................
382	c calcul de km et kn au debut du pas de temps
383	c.......................................................................
384	c
385	DO igrid=1,ngrid
386	kn(igrid,1)=knmin
387	km(igrid,1)=kmmin
388	kmpre(igrid,1)=km(igrid,1)
389	ENDDO
390	c
391	c-----------------------------------------------------------------------
392	DO ilev=2,nlev-1
393	DO igrid=1,ngrid
394	c-----------------------------------------------------------------------
395	c
396	kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
397	& *sn(igrid,ilev)
398	km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
399	& *sm(igrid,ilev)
400	kmpre(igrid,ilev)=km(igrid,ilev)
401	c
402	c-----------------------------------------------------------------------
403	ENDDO
404	ENDDO
405	c-----------------------------------------------------------------------
406	c
407	DO igrid=1,ngrid
408	kn(igrid,nlev)=kn(igrid,nlev-1)
409	km(igrid,nlev)=km(igrid,nlev-1)
410	kmpre(igrid,nlev)=km(igrid,nlev)
411	ENDDO
412	c
413	c.......................................................................
414	c boucle sur les niveaux 2 a nlev-1
415	c.......................................................................
416	c
417	c---->
418	DO 10001 ilev=2,nlev-1
419	c---->
420	DO 10002 igrid=1,ngrid
421	c
422	c.......................................................................
423	c
424	c calcul des termes sources et puits de l'equation de q2
425	c ------------------------------------------------------
426	c
427	knq3=kn(igrid,ilev)*snq2(igrid,ilev)
428	& /sn(igrid,ilev)
429	kmq3=km(igrid,ilev)*smq2(igrid,ilev)
430	& /sm(igrid,ilev)
431	c
432	termq=0.E+0
433	termq3=0.E+0
434	termqm2=0.E+0
435	termq3m2=0.E+0
436	c
437	tmp1=dt2.E+0 km(igrid,ilev)*m2(igrid,ilev)
438	tmp2=dt2.E+0 kmq3*m2(igrid,ilev)
439	termqm2=termqm2
440	& +dt2.E+0 km(igrid,ilev)*m2(igrid,ilev)
441	& -dt2.E+0 kmq3*m2(igrid,ilev)
442	termq3m2=termq3m2
443	& +dt2.E+0 kmq3*m2(igrid,ilev)
444	c
445	termq=termq
446	& -dt2.E+0 kn(igrid,ilev)*n2(igrid,ilev)
447	& +dt2.E+0 knq3*n2(igrid,ilev)
448	termq3=termq3
449	& -dt2.E+0 knq3*n2(igrid,ilev)
450	c
451	termq3=termq3
452	& -dt2.E+0 q(igrid,ilev)*3 / (b1long(igrid,ilev))
453	c
454	c.......................................................................
455	c
456	c resolution stationnaire couplee avec le gradient de vitesse local
457	c -----------------------------------------------------------------
458	c
459	c -----{on cherche le cisaillement qui annule l'equation de q^2
460	c supposee en q3}
461	c
462	tmp1=termq+termq3
463	tmp2=termqm2+termq3m2
464	m2cstat=m2(igrid,ilev)
465	& -(tmp1+tmp2)/(dt2.E+0km(igrid,ilev))
466	mcstat=sqrt(m2cstat)
467	c
468	c -----{puis on ecrit la valeur de q qui annule l'equation de m
469	c supposee en q3}
470	c
471	IF (ilev.eq.2) THEN
472	kmcstat=1.E+0 / mcstat
473	& ( unsdz(igrid,ilev)kmpre(igrid,ilev+1)
474	& *mpre(igrid,ilev+1)
475	& +unsdz(igrid,ilev-1)
476	& *cd(igrid)
477	& ( sqrt(u(igrid,3)2+v(igrid,3)*2)
478	& -mcstat/unsdzdec(igrid,ilev)
479	& -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
480	& /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
481	ELSE
482	kmcstat=1.E+0 / mcstat
483	& ( unsdz(igrid,ilev)kmpre(igrid,ilev+1)
484	& *mpre(igrid,ilev+1)
485	& +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
486	& *mpre(igrid,ilev-1) )
487	& /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
488	ENDIF
489	tmp2=kmcstat
490	& /( sm(igrid,ilev)/q2(igrid,ilev) )
491	& /long(igrid,ilev)
492	qcstat=tmp2**(1.E+0/3.E+0)
493	q2cstat=qcstat**2
494	c
495	c.......................................................................
496	c
497	c choix de la solution finale
498	c ---------------------------
499	c
500	q(igrid,ilev)=qcstat
501	q2(igrid,ilev)=q2cstat
502	m(igrid,ilev)=mcstat
503	m2(igrid,ilev)=m2cstat
504	c
505	c --->
506	c pour des raisons simples q2 est minore
507	c<---
508	c
509	IF (q2(igrid,ilev).lt.q2min) THEN
510	q2(igrid,ilev)=q2min
511	q(igrid,ilev)=sqrt(q2min)
512	ENDIF
513	c
514	c.......................................................................
515	c
516	c calcul final de kn et km
517	c ------------------------
518	c
519	gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
520	& * n2(igrid,ilev)
521	IF (gn.lt.gnmin) gn=gnmin
522	IF (gn.gt.gnmax) gn=gnmax
523	sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
524	sm(igrid,ilev)=
525	& (cm1+cm2*gn)
526	& /( (1.E+0 +cm3gn)(1.E+0 +cm4*gn) )
527	kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
528	& *sn(igrid,ilev)
529	km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
530	& *sm(igrid,ilev)
531	c
532	c.......................................................................
533	c
534	10002 CONTINUE
535	c
536	10001 CONTINUE
537	c
538	c.......................................................................
539	c
540	c
541	DO igrid=1,ngrid
542	kn(igrid,1)=knmin
543	km(igrid,1)=kmmin
544	q2(igrid,nlev)=q2(igrid,nlev-1)
545	q(igrid,nlev)=q(igrid,nlev-1)
546	kn(igrid,nlev)=kn(igrid,nlev-1)
547	km(igrid,nlev)=km(igrid,nlev-1)
548	ENDDO
549
550	c
551	c.......................................................................
552	c
553	RETURN
554	END

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

Download in other formats:

Original Format