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

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

Mars GCM: (a first step towards using parallel dynamics)

  • IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
    • advtrac.h (in dyn3d) removed and replaced by module infotrac.F
    • tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
  • Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
  • Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

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