Context Navigation

yamada4.F @ 1311

Last change on this file since 1311 was 1311, checked in by Laurent Fairhead, 14 years ago

Modifications to thermals for TKE transport

Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local

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

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

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

Download in other formats:

Original Format