Context Navigation

yamada4.F @ 1879

Last change on this file since 1879 was 1879, checked in by Ehouarn Millour, 11 years ago
Added some missing initializations. EM+JG
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 20.8 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
5	s ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
6	s ,iflag_pbl)
7	use dimphy
8	IMPLICIT NONE
9	#include "iniprint.h"
10	c.......................................................................
11	cym#include "dimensions.h"
12	cym#include "dimphy.h"
13	c.......................................................................
14	c
15	c dt : pas de temps
16	c g : g
17	c zlev : altitude a chaque niveau (interface inferieure de la couche
18	c de meme indice)
19	c zlay : altitude au centre de chaque couche
20	c u,v : vitesse au centre de chaque couche
21	c (en entree : la valeur au debut du pas de temps)
22	c teta : temperature potentielle au centre de chaque couche
23	c (en entree : la valeur au debut du pas de temps)
24	c cd : cdrag
25	c (en entree : la valeur au debut du pas de temps)
26	c q2 : $q^2$ au bas de chaque couche
27	c (en entree : la valeur au debut du pas de temps)
28	c (en sortie : la valeur a la fin du pas de temps)
29	c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
30	c couche)
31	c (en sortie : la valeur a la fin du pas de temps)
32	c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
33	c (en sortie : la valeur a la fin du pas de temps)
34	c
35	c iflag_pbl doit valoir entre 6 et 9
36	c l=6, on prend systematiquement une longueur d'equilibre
37	c iflag_pbl=6 : MY 2.0
38	c iflag_pbl=7 : MY 2.0.Fournier
39	c iflag_pbl=8/9 : MY 2.5
40	c iflag_pbl=8 with special obsolete treatments for convergence
41	c with Cmpi5 NPv3.1 simulations
42	c iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact
43	c iflag_pbl=12 = 11 with vertical diffusion off q2
44	c
45	c 2013/04/01 (FH hourdin@lmd.jussieu.fr)
46	c Correction for very stable PBLs (iflag_pbl=10 and 11)
47	c iflag_pbl=8 converges numerically with NPv3.1
48	c iflag_pbl=11 -> the model starts with NP from start files created by ce0l
49	c -> the model can run with longer time-steps.
50	c.......................................................................
51
52	REAL dt,g,rconst
53	real plev(klon,klev+1),temp(klon,klev)
54	real ustar(klon)
55	real kmin,qmin,pblhmin(klon),coriol(klon)
56	REAL zlev(klon,klev+1)
57	REAL zlay(klon,klev)
58	REAL u(klon,klev)
59	REAL v(klon,klev)
60	REAL teta(klon,klev)
61	REAL cd(klon)
62	REAL q2(klon,klev+1),qpre
63	REAL unsdz(klon,klev)
64	REAL unsdzdec(klon,klev+1)
65
66	REAL km(klon,klev+1)
67	REAL kmpre(klon,klev+1),tmp2
68	REAL mpre(klon,klev+1)
69	REAL kn(klon,klev+1)
70	REAL kq(klon,klev+1)
71	real ff(klon,klev+1),delta(klon,klev+1)
72	real aa(klon,klev+1),aa0,aa1
73	integer iflag_pbl,ngrid
74	integer nlay,nlev
75
76	logical first
77	integer ipas
78	save first,ipas
79	cFH/IM data first,ipas/.true.,0/
80	data first,ipas/.false.,0/
81	c$OMP THREADPRIVATE( first,ipas)
82
83	integer ig,k
84
85
86	real ri,zrif,zalpha,zsm,zsn
87	real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
88
89	real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
90	real dtetadz(klon,klev+1)
91	real m2cstat,mcstat,kmcstat
92	real l(klon,klev+1)
93	real,allocatable,save :: l0(:)
94	c$OMP THREADPRIVATE(l0)
95	real sq(klon),sqz(klon),zz(klon,klev+1)
96	integer iter
97
98	real ric,rifc,b1,kap
99	save ric,rifc,b1,kap
100	data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
101	c$OMP THREADPRIVATE(ric,rifc,b1,kap)
102	real frif,falpha,fsm
103	real fl,zzz,zl0,zq2,zn2
104
105	real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
106	s ,lyam(klon,klev),knyam(klon,klev)
107	s ,w2yam(klon,klev),t2yam(klon,klev)
108	logical,save :: firstcall=.true.
109	c$OMP THREADPRIVATE(firstcall)
110	frif(ri)=0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
111	falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
112	fsm(ri)=1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
113	fl(zzz,zl0,zq2,zn2)=
114	s max(min(l0(ig)kapzlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
115	s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
116
117
118	nlay=klev
119	nlev=klev+1
120
121	if (firstcall) then
122	allocate(l0(klon))
123	firstcall=.false.
124	endif
125
126
127	if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.12)) then
128	stop'probleme de coherence dans appel a MY'
129	endif
130
131	ipas=ipas+1
132
133
134	c.......................................................................
135	c les increments verticaux
136	c.......................................................................
137	c
138	c!!!!! allerte !!!!!c
139	c!!!!! zlev n'est pas declare a nlev !!!!!c
140	c!!!!! ---->
141	DO ig=1,ngrid
142	zlev(ig,nlev)=zlay(ig,nlay)
143	& +( zlay(ig,nlay) - zlev(ig,nlev-1) )
144	ENDDO
145	c!!!!! <----
146	c!!!!! allerte !!!!!c
147	c
148	DO k=1,nlay
149	DO ig=1,ngrid
150	unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
151	ENDDO
152	ENDDO
153	DO ig=1,ngrid
154	unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
155	ENDDO
156	DO k=2,nlay
157	DO ig=1,ngrid
158	unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
159	ENDDO
160	ENDDO
161	DO ig=1,ngrid
162	unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
163	ENDDO
164	c
165	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
166	! Computing M^2, N^2, Richardson numbers, stability functions
167	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
168	! initialize arrays:
169	m2(:,:)=0.0
170	sm(:,:)=0.0
171	rif(:,:)=0.0
172
173	do k=2,klev
174	do ig=1,ngrid
175	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
176	m2(ig,k)=((u(ig,k)-u(ig,k-1))2+(v(ig,k)-v(ig,k-1))2)
177	s /(dz(ig,k)*dz(ig,k))
178	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
179	n2(ig,k)=g2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
180	c n2(ig,k)=0.
181	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
182	if (ri.lt.ric) then
183	rif(ig,k)=frif(ri)
184	else
185	rif(ig,k)=rifc
186	endif
187	if(rif(ig,k).lt.0.16) then
188	alpha(ig,k)=falpha(rif(ig,k))
189	sm(ig,k)=fsm(rif(ig,k))
190	else
191	alpha(ig,k)=1.12
192	sm(ig,k)=0.085
193	endif
194	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
195	enddo
196	enddo
197
198
199	c====================================================================
200	c Computing the mixing length
201	c====================================================================
202
203	c Mise a jour de l0
204	if (iflag_pbl==8.or.iflag_pbl==10) then
205
206	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
207	! Iterative computation of l0
208	! This version is kept for iflag_pbl only for convergence
209	! with NPv3.1 Cmip5 simulations
210	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
211
212	do ig=1,ngrid
213	sq(ig)=1.e-10
214	sqz(ig)=1.e-10
215	enddo
216	do k=2,klev-1
217	do ig=1,ngrid
218	zq=sqrt(q2(ig,k))
219	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
220	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
221	enddo
222	enddo
223	do ig=1,ngrid
224	l0(ig)=0.2*sqz(ig)/sq(ig)
225	enddo
226	do k=2,klev
227	do ig=1,ngrid
228	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
229	enddo
230	enddo
231	! print*,'L0 cas 8 ou 10 ',l0
232
233	else
234
235	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
236	! In all other case, the assymptotic mixing length l0 is imposed (100m)
237	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
238
239	l0(:)=150.
240	do k=2,klev
241	do ig=1,ngrid
242	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
243	enddo
244	enddo
245	! print*,'L0 cas autres ',l0
246
247	endif
248
249
250	c====================================================================
251	c Yamada 2.0
252	c====================================================================
253	if (iflag_pbl.eq.6) then
254
255	do k=2,klev
256	q2(:,k)=l(:,k)*2zz(:,k)
257	enddo
258
259
260	else if (iflag_pbl.eq.7) then
261	c====================================================================
262	c Yamada 2.Fournier
263	c====================================================================
264
265	c Calcul de l, km, au pas precedent
266	do k=2,klev
267	do ig=1,ngrid
268	c print*,'SMML=',sm(ig,k),l(ig,k)
269	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
270	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
271	mpre(ig,k)=sqrt(m2(ig,k))
272	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
273	enddo
274	enddo
275
276	do k=2,klev-1
277	do ig=1,ngrid
278	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
279	mcstat=sqrt(m2cstat)
280
281	c print*,'M2 L=',k,mpre(ig,k),mcstat
282	c
283	c -----{puis on ecrit la valeur de q qui annule l'equation de m
284	c supposee en q3}
285	c
286	IF (k.eq.2) THEN
287	kmcstat=1.E+0 / mcstat
288	& ( unsdz(ig,k)kmpre(ig,k+1)
289	& *mpre(ig,k+1)
290	& +unsdz(ig,k-1)
291	& *cd(ig)
292	& ( sqrt(u(ig,3)2+v(ig,3)*2)
293	& -mcstat/unsdzdec(ig,k)
294	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
295	& /( unsdz(ig,k)+unsdz(ig,k-1) )
296	ELSE
297	kmcstat=1.E+0 / mcstat
298	& ( unsdz(ig,k)kmpre(ig,k+1)
299	& *mpre(ig,k+1)
300	& +unsdz(ig,k-1)*kmpre(ig,k-1)
301	& *mpre(ig,k-1) )
302	& /( unsdz(ig,k)+unsdz(ig,k-1) )
303	ENDIF
304	c print*,'T2 L=',k,tmp2
305	tmp2=kmcstat
306	& /( sm(ig,k)/q2(ig,k) )
307	& /l(ig,k)
308	q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
309	c print*,'Q2 L=',k,q2(ig,k)
310	c
311	enddo
312	enddo
313
314	else if (iflag_pbl==8.or.iflag_pbl==9) then
315	c====================================================================
316	c Yamada 2.5 a la Didi
317	c====================================================================
318
319
320	c Calcul de l, km, au pas precedent
321	do k=2,klev
322	do ig=1,ngrid
323	c print*,'SMML=',sm(ig,k),l(ig,k)
324	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
325	if (delta(ig,k).lt.1.e-20) then
326	c print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
327	delta(ig,k)=1.e-20
328	endif
329	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
330	aa0=
331	s (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
332	aa1=
333	s (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
334	c abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
335	aa(ig,k)=aa1dt/(delta(ig,k)l(ig,k))
336	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
337	qpre=sqrt(q2(ig,k))
338	! if (iflag_pbl.eq.8 ) then
339	if (aa(ig,k).gt.0.) then
340	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
341	else
342	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
343	endif
344	! else ! iflag_pbl=9
345	! if (aa(ig,k)*qpre.gt.0.9) then
346	! q2(ig,k)=(qpre10.)*2
347	! else
348	! q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
349	! endif
350	! endif
351	q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
352	c print,'Q2 L=',k,q2(ig,k),qpreqpre
353	enddo
354	enddo
355
356	else if (iflag_pbl>=10) then
357
358	! print*,'Schema mixte D'
359	! print*,'Longueur ',l(:,:)
360	do k=2,klev-1
361	l(:,k)=max(l(:,k),1.)
362	km(:,k)=l(:,k)sqrt(q2(:,k))sm(:,k)
363	q2(:,k)=q2(:,k)+dtkm(:,k)m2(:,k)*(1.-rif(:,k))
364	q2(:,k)=min(max(q2(:,k),1.e-10),1.e4)
365	q2(:,k)=1./(1./sqrt(q2(:,k))+dt/(2l(:,k)b1))
366	q2(:,k)=q2(:,k)*q2(:,k)
367	enddo
368
369
370	else
371	stop'Cas nom prevu dans yamada4'
372
373	endif ! Fin du cas 8
374
375	c print*,'OK8'
376
377	c====================================================================
378	c Calcul des coefficients de m�ange
379	c====================================================================
380	do k=2,klev
381	c print*,'k=',k
382	do ig=1,ngrid
383	cabde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
384	zq=sqrt(q2(ig,k))
385	km(ig,k)=l(ig,k)zqsm(ig,k)
386	kn(ig,k)=km(ig,k)*alpha(ig,k)
387	kq(ig,k)=l(ig,k)zq0.2
388	c print*,'KML=',km(ig,k),kn(ig,k)
389	enddo
390	enddo
391	! initialize near-surface and top-layer mixing coefficients
392	kq(1:ngrid,1)=kq(1:ngrid,2) ! constant (ie no gradient) near the surface
393	kq(1:ngrid,klev+1)=0 ! zero at the top
394
395	! Transport diffusif vertical de la TKE.
396	if (iflag_pbl.ge.12) then
397	! print*,'YAMADA VDIF'
398	q2(:,1)=q2(:,2)
399	call vdif_q2(dt,g,rconst,ngrid,plev,temp,kq,q2)
400	endif
401
402	c Traitement des cas noctrunes avec l'introduction d'une longueur
403	c minilale.
404
405	c====================================================================
406	c Traitement particulier pour les cas tres stables.
407	c D'apres Holtslag Boville.
408
409	if (prt_level>1) THEN
410	print*,'YAMADA4 0'
411	endif !(prt_level>1) THEN
412	do ig=1,ngrid
413	coriol(ig)=1.e-4
414	pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
415	enddo
416
417	! print*,'pblhmin ',pblhmin
418	CTest a remettre 21 11 02
419	c test abd 13 05 02 if(0.eq.1) then
420	if(1==1) then
421	if(iflag_pbl==8.or.iflag_pbl==10) then
422
423	do k=2,klev
424	do ig=1,ngrid
425	if (teta(ig,2).gt.teta(ig,1)) then
426	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
427	kmin=kapzlev(ig,k)qmin
428	else
429	kmin=-1. ! kmin n'est utilise que pour les SL stables.
430	endif
431	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
432	c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
433	c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
434	kn(ig,k)=kmin
435	km(ig,k)=kmin
436	kq(ig,k)=kmin
437	c la longueur de melange est suposee etre l= kap z
438	c K=l q Sm d'ou q2=(K/l Sm)**2
439	q2(ig,k)=(qmin/sm(ig,k))**2
440	endif
441	enddo
442	enddo
443
444	else
445
446	do k=2,klev
447	do ig=1,ngrid
448	if (teta(ig,2).gt.teta(ig,1)) then
449	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
450	kmin=kapzlev(ig,k)qmin
451	else
452	kmin=-1. ! kmin n'est utilise que pour les SL stables.
453	endif
454	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
455	c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
456	c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
457	kn(ig,k)=kmin
458	km(ig,k)=kmin
459	kq(ig,k)=kmin
460	c la longueur de melange est suposee etre l= kap z
461	c K=l q Sm d'ou q2=(K/l Sm)**2
462	sm(ig,k)=1.
463	alpha(ig,k)=1.
464	q2(ig,k)=min((qmin/sm(ig,k))**2,10.)
465	zq=sqrt(q2(ig,k))
466	km(ig,k)=l(ig,k)zqsm(ig,k)
467	kn(ig,k)=km(ig,k)*alpha(ig,k)
468	kq(ig,k)=l(ig,k)zq0.2
469	endif
470	enddo
471	enddo
472	endif
473
474	endif
475
476	if (prt_level>1) THEN
477	print*,'YAMADA4 1'
478	endif !(prt_level>1) THEN
479	c Diagnostique pour stokage
480
481	if(1.eq.0)then
482	rino=rif
483	smyam(1:ngrid,1)=0.
484	styam(1:ngrid,1)=0.
485	lyam(1:ngrid,1)=0.
486	knyam(1:ngrid,1)=0.
487	w2yam(1:ngrid,1)=0.
488	t2yam(1:ngrid,1)=0.
489
490	smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
491	styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
492	lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
493	knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
494
495	c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
496
497	w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24
498	s +lyam(1:ngrid,2:klev)5.17kn(1:ngrid,2:klev)
499	s *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
500
501	t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev)
502	s dtetadz(1:ngrid,2:klev)*2
503	s /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
504	endif
505
506	c print*,'OKFIN'
507	first=.false.
508	return
509	end
510	SUBROUTINE vdif_q2(timestep,gravity,rconst,ngrid,plev,temp,
511	& kmy,q2)
512	use dimphy
513	IMPLICIT NONE
514	c.......................................................................
515	#include "dimensions.h"
516	cccc#include "dimphy.h"
517	c.......................................................................
518	c
519	c dt : pas de temps
520
521	real plev(klon,klev+1)
522	real temp(klon,klev)
523	real timestep
524	real gravity,rconst
525	real kstar(klon,klev+1),zz
526	real kmy(klon,klev+1)
527	real q2(klon,klev+1)
528	real deltap(klon,klev+1)
529	real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
530	integer ngrid
531
532	integer i,k
533
534	! print*,'RD=',rconst
535	do k=1,klev
536	do i=1,ngrid
537	c test
538	! print*,'i,k',i,k
539	! print*,'temp(i,k)=',temp(i,k)
540	! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
541	zz=(plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
542	kstar(i,k)=0.125(kmy(i,k+1)+kmy(i,k))zz*zz
543	s /(plev(i,k)-plev(i,k+1))*timestep
544	enddo
545	enddo
546
547	do k=2,klev
548	do i=1,ngrid
549	deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
550	enddo
551	enddo
552	do i=1,ngrid
553	deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
554	deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
555	denom(i,klev+1)=deltap(i,klev+1)+kstar(i,klev)
556	alpha(i,klev+1)=deltap(i,klev+1)*q2(i,klev+1)/denom(i,klev+1)
557	beta(i,klev+1)=kstar(i,klev)/denom(i,klev+1)
558	enddo
559
560	do k=klev,2,-1
561	do i=1,ngrid
562	denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))*
563	s kstar(i,k)+kstar(i,k-1)
564	c correction d'un bug 10 01 2001
565	alpha(i,k)=(q2(i,k)*deltap(i,k)
566	s +kstar(i,k)*alpha(i,k+1))/denom(i,k)
567	beta(i,k)=kstar(i,k-1)/denom(i,k)
568	enddo
569	enddo
570
571	c Si on recalcule q2(1)
572	if(1.eq.0) then
573	do i=1,ngrid
574	denom(i,1)=deltap(i,1)+(1-beta(i,2))*kstar(i,1)
575	q2(i,1)=(q2(i,1)*deltap(i,1)
576	s +kstar(i,1)*alpha(i,2))/denom(i,1)
577	enddo
578	endif
579	c sinon, on peut sauter cette boucle...
580
581	do k=2,klev+1
582	do i=1,ngrid
583	q2(i,k)=alpha(i,k)+beta(i,k)*q2(i,k-1)
584	enddo
585	enddo
586
587	return
588	end
589	SUBROUTINE vdif_q2e(timestep,gravity,rconst,ngrid,
590	& plev,temp,kmy,q2)
591	use dimphy
592	IMPLICIT NONE
593	c.......................................................................
594	#include "dimensions.h"
595	cccc#include "dimphy.h"
596	c.......................................................................
597	c
598	c dt : pas de temps
599
600	real plev(klon,klev+1)
601	real temp(klon,klev)
602	real timestep
603	real gravity,rconst
604	real kstar(klon,klev+1),zz
605	real kmy(klon,klev+1)
606	real q2(klon,klev+1)
607	real deltap(klon,klev+1)
608	real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
609	integer ngrid
610
611	integer i,k
612
613	do k=1,klev
614	do i=1,ngrid
615	zz=(plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
616	kstar(i,k)=0.125(kmy(i,k+1)+kmy(i,k))zz*zz
617	s /(plev(i,k)-plev(i,k+1))*timestep
618	enddo
619	enddo
620
621	do k=2,klev
622	do i=1,ngrid
623	deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
624	enddo
625	enddo
626	do i=1,ngrid
627	deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
628	deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
629	enddo
630
631	do k=klev,2,-1
632	do i=1,ngrid
633	q2(i,k)=q2(i,k)+
634	s ( kstar(i,k)*(q2(i,k+1)-q2(i,k))
635	s -kstar(i,k-1)*(q2(i,k)-q2(i,k-1)) )
636	s /deltap(i,k)
637	enddo
638	enddo
639
640	do i=1,ngrid
641	q2(i,1)=q2(i,1)+
642	s ( kstar(i,1)*(q2(i,2)-q2(i,1))
643	s )
644	s /deltap(i,1)
645	q2(i,klev+1)=q2(i,klev+1)+
646	s (
647	s -kstar(i,klev)*(q2(i,klev+1)-q2(i,klev)) )
648	s /deltap(i,klev+1)
649	enddo
650
651	return
652	end

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

Download in other formats:

Original Format