Context Navigation

vdif_kc.F @ 1047

Last change on this file since 1047 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

File size: 21.3 KB

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

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

Download in other formats:

Original Format