Context Navigation

yamada4.F @ 1738

Last change on this file since 1738 was 1738, checked in by idelkadi, 11 years ago

Modifications for numerical stability of the boundary layer parameterizations.
Concerning yamada4.F :

option with asymptotic mixing length l0 imposed (iflag_pbl=8 and 9)
option with a new temporal scheme (iflag_pbl=10 and 11)
Correction for very stable PBLs (iflag_pbl=10 and 11) iflag_pbl=8 converges numerically with NPv3.1 iflag_pbl=11 -> now starts with NP from start files created by ce0l

-> the model can run with longer time-steps.

Concerning thermals :

Introduction of an implicit computation of vertical advection in
the environment of thermal plumes in thermcell_dq
impl = 0 : explicit, 1 : implicit, -1 : old version
controled by iflag_thermals =

15, 16 run with impl=-1 : numerical convergence with NPv3
17, 18 run with impl=1 : more stable

15 and 17 correspond to the activation of the stratocumulus "bidouille"

Modified routines (phylmd/):
calltherm.F90 : for managing the various options of thermcell_dq
coef_diff_turb_mod.F90 : yamada4 called for iflag_pbl<= 18 instead of 11
physiq.F : desactivation of the vertical diffusion of TKE for iflag_pbl=10
thermcellV0_main.F90 : calling thermcell_dq with implicit scheme
thermcell_dq.F90 : thermcell_dq with optional explicit/implicit scheme
thermcell_main.F90 : various option for vertical transport by thermals
yamada4.F : Yamada scheme with new options

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 20.6 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
169	do k=2,klev
170	do ig=1,ngrid
171	dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
172	m2(ig,k)=((u(ig,k)-u(ig,k-1))2+(v(ig,k)-v(ig,k-1))2)
173	s /(dz(ig,k)*dz(ig,k))
174	dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
175	n2(ig,k)=g2.dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
176	c n2(ig,k)=0.
177	ri=n2(ig,k)/max(m2(ig,k),1.e-10)
178	if (ri.lt.ric) then
179	rif(ig,k)=frif(ri)
180	else
181	rif(ig,k)=rifc
182	endif
183	if(rif(ig,k).lt.0.16) then
184	alpha(ig,k)=falpha(rif(ig,k))
185	sm(ig,k)=fsm(rif(ig,k))
186	else
187	alpha(ig,k)=1.12
188	sm(ig,k)=0.085
189	endif
190	zz(ig,k)=b1m2(ig,k)(1.-rif(ig,k))*sm(ig,k)
191	enddo
192	enddo
193
194
195	c====================================================================
196	c Computing the mixing length
197	c====================================================================
198
199	c Mise a jour de l0
200	if (iflag_pbl==8.or.iflag_pbl==10) then
201
202	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
203	! Iterative computation of l0
204	! This version is kept for iflag_pbl only for convergence
205	! with NPv3.1 Cmip5 simulations
206	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
207
208	do ig=1,ngrid
209	sq(ig)=1.e-10
210	sqz(ig)=1.e-10
211	enddo
212	do k=2,klev-1
213	do ig=1,ngrid
214	zq=sqrt(q2(ig,k))
215	sqz(ig)=sqz(ig)+zqzlev(ig,k)(zlay(ig,k)-zlay(ig,k-1))
216	sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
217	enddo
218	enddo
219	do ig=1,ngrid
220	l0(ig)=0.2*sqz(ig)/sq(ig)
221	enddo
222	do k=2,klev
223	do ig=1,ngrid
224	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
225	enddo
226	enddo
227	! print*,'L0 cas 8 ou 10 ',l0
228
229	else
230
231	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
232	! In all other case, the assymptotic mixing length l0 is imposed (100m)
233	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
234
235	l0(:)=150.
236	do k=2,klev
237	do ig=1,ngrid
238	l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
239	enddo
240	enddo
241	! print*,'L0 cas autres ',l0
242
243	endif
244
245
246	c====================================================================
247	c Yamada 2.0
248	c====================================================================
249	if (iflag_pbl.eq.6) then
250
251	do k=2,klev
252	q2(:,k)=l(:,k)*2zz(:,k)
253	enddo
254
255
256	else if (iflag_pbl.eq.7) then
257	c====================================================================
258	c Yamada 2.Fournier
259	c====================================================================
260
261	c Calcul de l, km, au pas precedent
262	do k=2,klev
263	do ig=1,ngrid
264	c print*,'SMML=',sm(ig,k),l(ig,k)
265	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
266	kmpre(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
267	mpre(ig,k)=sqrt(m2(ig,k))
268	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
269	enddo
270	enddo
271
272	do k=2,klev-1
273	do ig=1,ngrid
274	m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
275	mcstat=sqrt(m2cstat)
276
277	c print*,'M2 L=',k,mpre(ig,k),mcstat
278	c
279	c -----{puis on ecrit la valeur de q qui annule l'equation de m
280	c supposee en q3}
281	c
282	IF (k.eq.2) THEN
283	kmcstat=1.E+0 / mcstat
284	& ( unsdz(ig,k)kmpre(ig,k+1)
285	& *mpre(ig,k+1)
286	& +unsdz(ig,k-1)
287	& *cd(ig)
288	& ( sqrt(u(ig,3)2+v(ig,3)*2)
289	& -mcstat/unsdzdec(ig,k)
290	& -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
291	& /( unsdz(ig,k)+unsdz(ig,k-1) )
292	ELSE
293	kmcstat=1.E+0 / mcstat
294	& ( unsdz(ig,k)kmpre(ig,k+1)
295	& *mpre(ig,k+1)
296	& +unsdz(ig,k-1)*kmpre(ig,k-1)
297	& *mpre(ig,k-1) )
298	& /( unsdz(ig,k)+unsdz(ig,k-1) )
299	ENDIF
300	c print*,'T2 L=',k,tmp2
301	tmp2=kmcstat
302	& /( sm(ig,k)/q2(ig,k) )
303	& /l(ig,k)
304	q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
305	c print*,'Q2 L=',k,q2(ig,k)
306	c
307	enddo
308	enddo
309
310	else if (iflag_pbl==8.or.iflag_pbl==9) then
311	c====================================================================
312	c Yamada 2.5 a la Didi
313	c====================================================================
314
315
316	c Calcul de l, km, au pas precedent
317	do k=2,klev
318	do ig=1,ngrid
319	c print*,'SMML=',sm(ig,k),l(ig,k)
320	delta(ig,k)=q2(ig,k)/(l(ig,k)*2sm(ig,k))
321	if (delta(ig,k).lt.1.e-20) then
322	c print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
323	delta(ig,k)=1.e-20
324	endif
325	km(ig,k)=l(ig,k)sqrt(q2(ig,k))sm(ig,k)
326	aa0=
327	s (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
328	aa1=
329	s (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
330	c abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
331	aa(ig,k)=aa1dt/(delta(ig,k)l(ig,k))
332	c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
333	qpre=sqrt(q2(ig,k))
334	! if (iflag_pbl.eq.8 ) then
335	if (aa(ig,k).gt.0.) then
336	q2(ig,k)=(qpre+aa(ig,k)qpreqpre)**2
337	else
338	q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
339	endif
340	! else ! iflag_pbl=9
341	! if (aa(ig,k)*qpre.gt.0.9) then
342	! q2(ig,k)=(qpre10.)*2
343	! else
344	! q2(ig,k)=(qpre/(1.-aa(ig,k)qpre))*2
345	! endif
346	! endif
347	q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
348	c print,'Q2 L=',k,q2(ig,k),qpreqpre
349	enddo
350	enddo
351
352	else if (iflag_pbl>=10) then
353
354	! print*,'Schema mixte D'
355	! print*,'Longueur ',l(:,:)
356	do k=2,klev-1
357	l(:,k)=max(l(:,k),1.)
358	km(:,k)=l(:,k)sqrt(q2(:,k))sm(:,k)
359	q2(:,k)=q2(:,k)+dtkm(:,k)m2(:,k)*(1.-rif(:,k))
360	q2(:,k)=min(max(q2(:,k),1.e-10),1.e4)
361	q2(:,k)=1./(1./sqrt(q2(:,k))+dt/(2l(:,k)b1))
362	q2(:,k)=q2(:,k)*q2(:,k)
363	enddo
364
365
366	else
367	stop'Cas nom prevu dans yamada4'
368
369	endif ! Fin du cas 8
370
371	c print*,'OK8'
372
373	c====================================================================
374	c Calcul des coefficients de m�ange
375	c====================================================================
376	do k=2,klev
377	c print*,'k=',k
378	do ig=1,ngrid
379	cabde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
380	zq=sqrt(q2(ig,k))
381	km(ig,k)=l(ig,k)zqsm(ig,k)
382	kn(ig,k)=km(ig,k)*alpha(ig,k)
383	kq(ig,k)=l(ig,k)zq0.2
384	c print*,'KML=',km(ig,k),kn(ig,k)
385	enddo
386	enddo
387
388	! Transport diffusif vertical de la TKE.
389	if (iflag_pbl.ge.12) then
390	! print*,'YAMADA VDIF'
391	q2(:,1)=q2(:,2)
392	call vdif_q2(dt,g,rconst,ngrid,plev,temp,kq,q2)
393	endif
394
395	c Traitement des cas noctrunes avec l'introduction d'une longueur
396	c minilale.
397
398	c====================================================================
399	c Traitement particulier pour les cas tres stables.
400	c D'apres Holtslag Boville.
401
402	if (prt_level>1) THEN
403	print*,'YAMADA4 0'
404	endif !(prt_level>1) THEN
405	do ig=1,ngrid
406	coriol(ig)=1.e-4
407	pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
408	enddo
409
410	! print*,'pblhmin ',pblhmin
411	CTest a remettre 21 11 02
412	c test abd 13 05 02 if(0.eq.1) then
413	if(1==1) then
414	if(iflag_pbl==8.or.iflag_pbl==10) then
415
416	do k=2,klev
417	do ig=1,ngrid
418	if (teta(ig,2).gt.teta(ig,1)) then
419	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
420	kmin=kapzlev(ig,k)qmin
421	else
422	kmin=-1. ! kmin n'est utilise que pour les SL stables.
423	endif
424	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
425	c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
426	c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
427	kn(ig,k)=kmin
428	km(ig,k)=kmin
429	kq(ig,k)=kmin
430	c la longueur de melange est suposee etre l= kap z
431	c K=l q Sm d'ou q2=(K/l Sm)**2
432	q2(ig,k)=(qmin/sm(ig,k))**2
433	endif
434	enddo
435	enddo
436
437	else
438
439	do k=2,klev
440	do ig=1,ngrid
441	if (teta(ig,2).gt.teta(ig,1)) then
442	qmin=ustar(ig)(max(1.-zlev(ig,k)/pblhmin(ig),0.))*2
443	kmin=kapzlev(ig,k)qmin
444	else
445	kmin=-1. ! kmin n'est utilise que pour les SL stables.
446	endif
447	if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
448	c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
449	c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
450	kn(ig,k)=kmin
451	km(ig,k)=kmin
452	kq(ig,k)=kmin
453	c la longueur de melange est suposee etre l= kap z
454	c K=l q Sm d'ou q2=(K/l Sm)**2
455	sm(ig,k)=1.
456	alpha(ig,k)=1.
457	q2(ig,k)=min((qmin/sm(ig,k))**2,10.)
458	zq=sqrt(q2(ig,k))
459	km(ig,k)=l(ig,k)zqsm(ig,k)
460	kn(ig,k)=km(ig,k)*alpha(ig,k)
461	kq(ig,k)=l(ig,k)zq0.2
462	endif
463	enddo
464	enddo
465	endif
466
467	endif
468
469	if (prt_level>1) THEN
470	print*,'YAMADA4 1'
471	endif !(prt_level>1) THEN
472	c Diagnostique pour stokage
473
474	if(1.eq.0)then
475	rino=rif
476	smyam(1:ngrid,1)=0.
477	styam(1:ngrid,1)=0.
478	lyam(1:ngrid,1)=0.
479	knyam(1:ngrid,1)=0.
480	w2yam(1:ngrid,1)=0.
481	t2yam(1:ngrid,1)=0.
482
483	smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
484	styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
485	lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
486	knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
487
488	c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
489
490	w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24
491	s +lyam(1:ngrid,2:klev)5.17kn(1:ngrid,2:klev)
492	s *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
493
494	t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev)
495	s dtetadz(1:ngrid,2:klev)*2
496	s /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
497	endif
498
499	c print*,'OKFIN'
500	first=.false.
501	return
502	end
503	SUBROUTINE vdif_q2(timestep,gravity,rconst,ngrid,plev,temp,
504	& kmy,q2)
505	use dimphy
506	IMPLICIT NONE
507	c.......................................................................
508	#include "dimensions.h"
509	cccc#include "dimphy.h"
510	c.......................................................................
511	c
512	c dt : pas de temps
513
514	real plev(klon,klev+1)
515	real temp(klon,klev)
516	real timestep
517	real gravity,rconst
518	real kstar(klon,klev+1),zz
519	real kmy(klon,klev+1)
520	real q2(klon,klev+1)
521	real deltap(klon,klev+1)
522	real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
523	integer ngrid
524
525	integer i,k
526
527	! print*,'RD=',rconst
528	do k=1,klev
529	do i=1,ngrid
530	c test
531	! print*,'i,k',i,k
532	! print*,'temp(i,k)=',temp(i,k)
533	! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
534	zz=(plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
535	kstar(i,k)=0.125(kmy(i,k+1)+kmy(i,k))zz*zz
536	s /(plev(i,k)-plev(i,k+1))*timestep
537	enddo
538	enddo
539
540	do k=2,klev
541	do i=1,ngrid
542	deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
543	enddo
544	enddo
545	do i=1,ngrid
546	deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
547	deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
548	denom(i,klev+1)=deltap(i,klev+1)+kstar(i,klev)
549	alpha(i,klev+1)=deltap(i,klev+1)*q2(i,klev+1)/denom(i,klev+1)
550	beta(i,klev+1)=kstar(i,klev)/denom(i,klev+1)
551	enddo
552
553	do k=klev,2,-1
554	do i=1,ngrid
555	denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))*
556	s kstar(i,k)+kstar(i,k-1)
557	c correction d'un bug 10 01 2001
558	alpha(i,k)=(q2(i,k)*deltap(i,k)
559	s +kstar(i,k)*alpha(i,k+1))/denom(i,k)
560	beta(i,k)=kstar(i,k-1)/denom(i,k)
561	enddo
562	enddo
563
564	c Si on recalcule q2(1)
565	if(1.eq.0) then
566	do i=1,ngrid
567	denom(i,1)=deltap(i,1)+(1-beta(i,2))*kstar(i,1)
568	q2(i,1)=(q2(i,1)*deltap(i,1)
569	s +kstar(i,1)*alpha(i,2))/denom(i,1)
570	enddo
571	endif
572	c sinon, on peut sauter cette boucle...
573
574	do k=2,klev+1
575	do i=1,ngrid
576	q2(i,k)=alpha(i,k)+beta(i,k)*q2(i,k-1)
577	enddo
578	enddo
579
580	return
581	end
582	SUBROUTINE vdif_q2e(timestep,gravity,rconst,ngrid,
583	& plev,temp,kmy,q2)
584	use dimphy
585	IMPLICIT NONE
586	c.......................................................................
587	#include "dimensions.h"
588	cccc#include "dimphy.h"
589	c.......................................................................
590	c
591	c dt : pas de temps
592
593	real plev(klon,klev+1)
594	real temp(klon,klev)
595	real timestep
596	real gravity,rconst
597	real kstar(klon,klev+1),zz
598	real kmy(klon,klev+1)
599	real q2(klon,klev+1)
600	real deltap(klon,klev+1)
601	real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
602	integer ngrid
603
604	integer i,k
605
606	do k=1,klev
607	do i=1,ngrid
608	zz=(plev(i,k)+plev(i,k+1))gravity/(rconsttemp(i,k))
609	kstar(i,k)=0.125(kmy(i,k+1)+kmy(i,k))zz*zz
610	s /(plev(i,k)-plev(i,k+1))*timestep
611	enddo
612	enddo
613
614	do k=2,klev
615	do i=1,ngrid
616	deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
617	enddo
618	enddo
619	do i=1,ngrid
620	deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
621	deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
622	enddo
623
624	do k=klev,2,-1
625	do i=1,ngrid
626	q2(i,k)=q2(i,k)+
627	s ( kstar(i,k)*(q2(i,k+1)-q2(i,k))
628	s -kstar(i,k-1)*(q2(i,k)-q2(i,k-1)) )
629	s /deltap(i,k)
630	enddo
631	enddo
632
633	do i=1,ngrid
634	q2(i,1)=q2(i,1)+
635	s ( kstar(i,1)*(q2(i,2)-q2(i,1))
636	s )
637	s /deltap(i,1)
638	q2(i,klev+1)=q2(i,klev+1)+
639	s (
640	s -kstar(i,klev)*(q2(i,klev+1)-q2(i,klev)) )
641	s /deltap(i,klev+1)
642	enddo
643
644	return
645	end

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

Download in other formats:

Original Format