source: trunk/LMDZ.MARS/libf/phymars/vdif_kc.F @ 706

Last change on this file since 706 was 325, checked in by acolaitis, 13 years ago

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