Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

vdif_kcay.f90 @ 5354

Last change on this file since 5354 was 5268, checked in by abarral, 7 weeks ago
.f90 <-> .F90 depending on cpp key use
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.3 KB

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

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

Download in other formats:

Original Format