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

Last change on this file since 3026 was 2823, checked in by emillour, 2 years ago

Mars GCM:
Remove the "tracer" (logical) flag as we now always run with at least
one tracer.
EM
File size: 22.4 KB
RevLine 
[1047]1      SUBROUTINE vdif_kc(ngrid,nlay,nq,dt,g,
2     &                   zlev,zlay,u,v,teta,cd,q2,km,kn,zq)
[1036]3
[1047]4      use tracer_mod, only: noms
[38]5      IMPLICIT NONE
[2823]6
[38]7c.......................................................................
8c
9c dt : pas de temps
10c g  : g
11c zlev : altitude a chaque niveau (interface inferieure de la couche
12c        de meme indice)
13c zlay : altitude au centre de chaque couche
14c u,v : vitesse au centre de chaque couche
15c       (en entree : la valeur au debut du pas de temps)
16c teta : temperature potentielle au centre de chaque couche
17c        (en entree : la valeur au debut du pas de temps)
18c cd : cdrag
19c      (en entree : la valeur au debut du pas de temps)
20c q2 : $q^2$ au bas de chaque couche
21c      (en entree : la valeur au debut du pas de temps)
22c      (en sortie : la valeur a la fin du pas de temps)
23c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
24c      couche)
25c      (en sortie : la valeur a la fin du pas de temps)
26c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
27c      (en sortie : la valeur a la fin du pas de temps)
28c
29c.......................................................................
[1130]30      REAL,INTENT(IN) :: dt,g
31      REAL,INTENT(IN) :: zlev(ngrid,nlay+1)
32      REAL,INTENT(IN) :: zlay(ngrid,nlay)
33      REAL,INTENT(IN) :: u(ngrid,nlay)
34      REAL,INTENT(IN) :: v(ngrid,nlay)
35      REAL,INTENT(IN) :: teta(ngrid,nlay)
36      REAL,INTENT(IN) :: cd(ngrid)
37      REAL,INTENT(INOUT) :: q2(ngrid,nlay+1)
38      REAL,INTENT(OUT) :: km(ngrid,nlay+1)
39      REAL,INTENT(OUT) :: kn(ngrid,nlay+1)
40      REAL,INTENT(IN) :: zq(ngrid,nlay,nq)
[38]41c.......................................................................
42c
43c nlay : nombre de couches       
44c nlev : nombre de niveaux
45c ngrid : nombre de points de grille       
46c unsdz : 1 sur l'epaisseur de couche
47c unsdzdec : 1 sur la distance entre le centre de la couche et le
48c            centre de la couche inferieure
49c q : echelle de vitesse au bas de chaque couche
50c     (valeur a la fin du pas de temps)
51c
52c.......................................................................
[1047]53      INTEGER,INTENT(IN) :: nlay,ngrid,nq
54      INTEGER nlev
55      REAL unsdz(ngrid,nlay)
56      REAL unsdzdec(ngrid,nlay+1)
57      REAL q(ngrid,nlay+1)
[38]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.......................................................................
[1047]67      REAL kmpre(ngrid,nlay+1)
[38]68      REAL qcstat
69      REAL q2cstat
70c.......................................................................
71c
72c long : longueur de melange calculee selon Blackadar
73c
74c.......................................................................
[1047]75      REAL long(ngrid,nlay+1)
[38]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
[1047]98      REAL m(ngrid,nlay+1)
99      REAL mpre(ngrid,nlay+1)
100      REAL m2(ngrid,nlay+1)
101      REAL n2(ngrid,nlay+1)
[38]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
[1047]123c      REAL ri(ngrid,nlay+1)
124      REAL sn(ngrid,nlay+1)
125      REAL snq2(ngrid,nlay+1)
126      REAL sm(ngrid,nlay+1)
127      REAL smq2(ngrid,nlay+1)
[38]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)
186c
187      PARAMETER (
188     &  cn1=a2*(1.E+0 -6.E+0 *a1/b1)
189     &          )
190      PARAMETER (
191     &  cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
192     &          )
193      PARAMETER (
194     &  cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
195     &          )
196      PARAMETER (
197     &  cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*(1.E+0 -6.E+0 *a1/b1)
198     &          -3.E+0 *c1*(b2+6.E+0 *a1)))
199     &          )
200      PARAMETER (
201     &  cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
202     &          )
203      PARAMETER (
204     &  cm4=-9.E+0 *a1*a2
205     &          )
[325]206
207c AC: variables for theta_m computation
208
209      INTEGER ico2,iq
210      SAVE ico2
[2616]211
212!$OMP THREADPRIVATE(ico2)
213
[325]214      REAL m_co2, m_noco2, A , B
215      SAVE A, B
[2616]216
217!$OMP THREADPRIVATE(A,B)
218
[325]219      LOGICAL firstcall
220      save firstcall
[2616]221
222!$OMP THREADPRIVATE(firstcall)
223
[325]224      data firstcall/.true./
[1047]225      REAL zhc(ngrid,nlay)
[38]226c.......................................................................
[325]227c  Initialization
228c.......................................................................
229
[1779]230      ! AS: OK firstcall absolute
[325]231      if(firstcall) then
232        ico2=0
233!     Prepare Special treatment if one of the tracers is CO2 gas
[1047]234           do iq=1,nq
[325]235             if (noms(iq).eq."co2") then
236                ico2=iq
237                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
238                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
239!               Compute A and B coefficient use to compute
240!               mean molecular mass Mair defined by
241!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
242!               1/Mair = A*q(ico2) + B
243                A =(1/m_co2 - 1/m_noco2)
244                B=1/m_noco2
245             end if
246           enddo
247
248      firstcall=.false.
249      endif !of if firstcall
250
[1047]251      nlev=nlay+1
252     
[1130]253      ! Other initializations of local variables (to be clean)
254      long(:,:)=0.
255      n2(:,:)=0.
256      sn(:,:)=0.
257      snq2(:,:)=0.
258      sm(:,:)=0.
259      smq2(:,:)=0.
[325]260c.......................................................................
261c  Special treatment for co2
262c.......................................................................
263
264      if (ico2.ne.0) then
265!     Special case if one of the tracers is CO2 gas
266         DO ilay=1,nlay
267           DO igrid=1,ngrid
268            zhc(igrid,ilay) = teta(igrid,ilay)*(A*zq(igrid,ilay,ico2)+B)
269           ENDDO
270         ENDDO
271       else
272          zhc(:,:)=teta(:,:)
273       end if
274
275c.......................................................................
[38]276c  traitment des valeur de q2 en entree
277c.......................................................................
278c
279      DO ilev=1,nlev
280                                                      DO igrid=1,ngrid
281        q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
282        q(igrid,ilev)=sqrt(q2(igrid,ilev))
283                                                      ENDDO
284      ENDDO
285c
286                                                      DO igrid=1,ngrid
287      tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
288      q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
289      q2(igrid,1)=amax1(q2(igrid,1),q2min)
290      q(igrid,1)=sqrt(q2(igrid,1))
291                                                      ENDDO
292c
293c.......................................................................
294c  les increments verticaux
295c.......................................................................
296c
297c!!!!! allerte !!!!!c
298c!!!!! zlev n'est pas declare a nlev !!!!!c
299c!!!!! ---->
300c                                                     DO igrid=1,ngrid
301c           zlev(igrid,nlev)=zlay(igrid,nlay)
302c    &             +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
303c                                                     ENDDO           
304c!!!!! <----
305c!!!!! allerte !!!!!c
306c
307      DO ilay=1,nlay
308                                                      DO igrid=1,ngrid
309        unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
310                                                      ENDDO
311      ENDDO
312                                                      DO igrid=1,ngrid
313      unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
314                                                      ENDDO
315      DO ilay=2,nlay
316                                                      DO igrid=1,ngrid
317        unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
318                                                      ENDDO
319      ENDDO
320                                                      DO igrid=1,ngrid
321      unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
322                                                      ENDDO
323c
324c.......................................................................
325c  le cisaillement et le gradient de temperature
326c.......................................................................
327c
328                                                      DO igrid=1,ngrid
329      m2(igrid,1)=(unsdzdec(igrid,1)
330     &                   *u(igrid,1))**2
331     &                 +(unsdzdec(igrid,1)
332     &                   *v(igrid,1))**2
333      m(igrid,1)=sqrt(m2(igrid,1))
334      mpre(igrid,1)=m(igrid,1)
335                                                      ENDDO
336c
337c-----------------------------------------------------------------------
338      DO ilev=2,nlev-1
339                                                      DO igrid=1,ngrid
340c-----------------------------------------------------------------------
341c
342        n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
[325]343     &                   *(zhc(igrid,ilev)-zhc(igrid,ilev-1))
344     &                   /(zhc(igrid,ilev)+zhc(igrid,ilev-1)) *2.E+0
[38]345c
346c --->
347c       on ne sais traiter que les cas stratifies. et l'ajustement
348c       convectif est cense faire en sorte que seul des configurations
349c       stratifiees soient rencontrees en entree de cette routine.
350c       mais, bon ... on sait jamais (meme on sait que n2 prends
351c       quelques valeurs negatives ... parfois) alors :
352c<---
353c
354        IF (n2(igrid,ilev).lt.0.E+0) THEN
355          n2(igrid,ilev)=0.E+0
356        ENDIF
357c
358        m2(igrid,ilev)=(unsdzdec(igrid,ilev)
359     &                     *(u(igrid,ilev)-u(igrid,ilev-1)))**2
360     &                   +(unsdzdec(igrid,ilev)
361     &                     *(v(igrid,ilev)-v(igrid,ilev-1)))**2
362        m(igrid,ilev)=sqrt(m2(igrid,ilev))
363        mpre(igrid,ilev)=m(igrid,ilev)
364c
365c-----------------------------------------------------------------------
366                                                      ENDDO
367      ENDDO
368c-----------------------------------------------------------------------
369c
370                                                      DO igrid=1,ngrid
371      m2(igrid,nlev)=m2(igrid,nlev-1)
372      m(igrid,nlev)=m(igrid,nlev-1)
373      mpre(igrid,nlev)=m(igrid,nlev)
374                                                      ENDDO
375c
376c.......................................................................
377c  calcul des fonctions de stabilite
378c.......................................................................
379c
380c-----------------------------------------------------------------------
381      DO ilev=2,nlev-1
382                                                      DO igrid=1,ngrid
383c-----------------------------------------------------------------------
384c
385        tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
386        long(igrid,ilev)=tmp1/(1.E+0 + tmp1/long0)
387        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
388     &                                           * n2(igrid,ilev)
389        gm=long(igrid,ilev)**2 / q2(igrid,ilev)
390     &                                           * m2(igrid,ilev)
391c
392        gninf=.false.
393        gnsup=.false.
[1130]394        long(igrid,ilev)=long(igrid,ilev) ! not very useful...
395        long(igrid,ilev)=long(igrid,ilev) ! not very useful...
[38]396c
397        IF (gn.lt.gnmin) THEN
398          gninf=.true.
399          gn=gnmin
400        ENDIF
401c
402        IF (gn.gt.gnmax) THEN
403          gnsup=.true.
404          gn=gnmax
405        ENDIF
406c
407        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
408        sm(igrid,ilev)=
409     &    (cm1+cm2*gn)
410     &   /( (1.E+0 +cm3*gn)
411     &     *(1.E+0 +cm4*gn) )
412c
413        IF ((gninf).or.(gnsup)) THEN
414          snq2(igrid,ilev)=0.E+0
415          smq2(igrid,ilev)=0.E+0
416        ELSE
417          snq2(igrid,ilev)=
418     &     -gn
419     &     *(-cn1*cn2/(1.E+0 +cn2*gn)**2 )
420          smq2(igrid,ilev)=
421     &     -gn
422     &     *( cm2*(1.E+0 +cm3*gn)
423     &           *(1.E+0 +cm4*gn)
424     &       -( cm3*(1.E+0 +cm4*gn)
425     &         +cm4*(1.E+0 +cm3*gn) )
426     &       *(cm1+cm2*gn)            )
427     &     /( (1.E+0 +cm3*gn)
428     &       *(1.E+0 +cm4*gn) )**2
429        ENDIF
430c
431c --->
432c       la decomposition de Taylor en q2 n'a de sens que
433c       dans les cas stratifies ou sn et sm sont quasi
434c       proportionnels a q2. ailleurs on laisse le meme
435c       algorithme car l'ajustement convectif fait le travail.
436c       mais c'est delirant quand sn et snq2 n'ont pas le meme
437c       signe : dans ces cas, on ne fait pas la decomposition.
438c<---
439c
440        IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.E+0)
441     &      snq2(igrid,ilev)=0.E+0
442        IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.E+0)
443     &      smq2(igrid,ilev)=0.E+0
444c
445c-----------------------------------------------------------------------
446                                                      ENDDO
447      ENDDO
448c-----------------------------------------------------------------------
449c
450c.......................................................................
451c  calcul de km et kn au debut du pas de temps
452c.......................................................................
453c
454                                                      DO igrid=1,ngrid
455      kn(igrid,1)=knmin
456      km(igrid,1)=kmmin
457      kmpre(igrid,1)=km(igrid,1)
458                                                      ENDDO
459c
460c-----------------------------------------------------------------------
461      DO ilev=2,nlev-1
462                                                      DO igrid=1,ngrid
463c-----------------------------------------------------------------------
464c
465        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
466     &                                         *sn(igrid,ilev)
467        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
468     &                                         *sm(igrid,ilev)
469        kmpre(igrid,ilev)=km(igrid,ilev)
470c
471c-----------------------------------------------------------------------
472                                                      ENDDO
473      ENDDO
474c-----------------------------------------------------------------------
475c
476                                                      DO igrid=1,ngrid
477      kn(igrid,nlev)=kn(igrid,nlev-1)
478      km(igrid,nlev)=km(igrid,nlev-1)
479      kmpre(igrid,nlev)=km(igrid,nlev)
480                                                      ENDDO
481c
482c.......................................................................
483c  boucle sur les niveaux 2 a nlev-1
484c.......................................................................
485c
486c---->
487      DO 10001 ilev=2,nlev-1
488c---->
489      DO 10002 igrid=1,ngrid
490c
491c.......................................................................
492c
493c  calcul des termes sources et puits de l'equation de q2
494c  ------------------------------------------------------
495c
496        knq3=kn(igrid,ilev)*snq2(igrid,ilev)
497     &                                    /sn(igrid,ilev)
498        kmq3=km(igrid,ilev)*smq2(igrid,ilev)
499     &                                    /sm(igrid,ilev)
500c
501        termq=0.E+0
502        termq3=0.E+0
503        termqm2=0.E+0
504        termq3m2=0.E+0
505c
506        tmp1=dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
507        tmp2=dt*2.E+0 *kmq3*m2(igrid,ilev)
508        termqm2=termqm2
509     &    +dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
510     &    -dt*2.E+0 *kmq3*m2(igrid,ilev)
511        termq3m2=termq3m2
512     &    +dt*2.E+0 *kmq3*m2(igrid,ilev)
513c
514        termq=termq
515     &    -dt*2.E+0 *kn(igrid,ilev)*n2(igrid,ilev)
516     &    +dt*2.E+0 *knq3*n2(igrid,ilev)
517        termq3=termq3
518     &    -dt*2.E+0 *knq3*n2(igrid,ilev)
519c
520        termq3=termq3
521     &    -dt*2.E+0 *q(igrid,ilev)**3 / (b1*long(igrid,ilev))
522c
523c.......................................................................
524c
525c  resolution stationnaire couplee avec le gradient de vitesse local
526c  -----------------------------------------------------------------
527c
528c  -----{on cherche le cisaillement qui annule l'equation de q^2
529c        supposee en q3}
530c
531        tmp1=termq+termq3
532        tmp2=termqm2+termq3m2
533        m2cstat=m2(igrid,ilev)
534     &      -(tmp1+tmp2)/(dt*2.E+0*km(igrid,ilev))
535        mcstat=sqrt(m2cstat)
536c
537c  -----{puis on ecrit la valeur de q qui annule l'equation de m
538c        supposee en q3}
539c
540        IF (ilev.eq.2) THEN
541          kmcstat=1.E+0 / mcstat
542     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
543     &                        *mpre(igrid,ilev+1)
544     &      +unsdz(igrid,ilev-1)
545     &              *cd(igrid)
546     &              *( sqrt(u(igrid,3)**2+v(igrid,3)**2)
547     &                -mcstat/unsdzdec(igrid,ilev)
548     &                -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
549     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
550        ELSE
551          kmcstat=1.E+0 / mcstat
552     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
553     &                        *mpre(igrid,ilev+1)
554     &      +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
555     &                          *mpre(igrid,ilev-1) )
556     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
557        ENDIF
558        tmp2=kmcstat
559     &      /( sm(igrid,ilev)/q2(igrid,ilev) )
560     &      /long(igrid,ilev)
561        qcstat=tmp2**(1.E+0/3.E+0)
562        q2cstat=qcstat**2
563c
564c.......................................................................
565c
566c  choix de la solution finale
567c  ---------------------------
568c
569          q(igrid,ilev)=qcstat
570          q2(igrid,ilev)=q2cstat
571          m(igrid,ilev)=mcstat
572          m2(igrid,ilev)=m2cstat
573c
574c --->
575c       pour des raisons simples q2 est minore
576c<---
577c
578        IF (q2(igrid,ilev).lt.q2min) THEN
579          q2(igrid,ilev)=q2min
580          q(igrid,ilev)=sqrt(q2min)
581        ENDIF
582c
583c.......................................................................
584c
585c  calcul final de kn et km
586c  ------------------------
587c
588        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
589     &                                           * n2(igrid,ilev)
590        IF (gn.lt.gnmin) gn=gnmin
591        IF (gn.gt.gnmax) gn=gnmax
592        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
593        sm(igrid,ilev)=
594     &    (cm1+cm2*gn)
595     &   /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
596        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
597     &                 *sn(igrid,ilev)
598        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
599     &                 *sm(igrid,ilev)
600c
601c.......................................................................
602c
60310002 CONTINUE
604c
60510001 CONTINUE
606c
607c.......................................................................
608c
609c
610                                                      DO igrid=1,ngrid
611      kn(igrid,1)=knmin
612      km(igrid,1)=kmmin
613      q2(igrid,nlev)=q2(igrid,nlev-1)
614      q(igrid,nlev)=q(igrid,nlev-1)
615      kn(igrid,nlev)=kn(igrid,nlev-1)
616      km(igrid,nlev)=km(igrid,nlev-1)
617                                                      ENDDO
618
619c
620c.......................................................................
621c
[1130]622
623!        call writediagfi(ngrid,'vdif_kc_q2','','',3,q2(:,1:nlay))
624!        call writediagfi(ngrid,'vdif_kc_km','','',3,km(:,1:nlay))
625!        call writediagfi(ngrid,'vdif_kc_kn','','',3,kn(:,1:nlay))
626!        call writediagfi(ngrid,'vdif_kc_unsdz','','',3,unsdz(:,1:nlay))
627!        call writediagfi(ngrid,'vdif_kc_unsddecz','','',3,
628!     &                   unsdzdec(:,1:nlay))
629!        call writediagfi(ngrid,'vdif_kc_q','','',3,q(:,1:nlay))
630!        call writediagfi(ngrid,'vdif_kc_kmpre','','',3,
631!     &                   kmpre(:,1:nlay))
632!        call writediagfi(ngrid,'vdif_kc_long','','',3,long(:,1:nlay))
633!        call writediagfi(ngrid,'vdif_kc_sn','','',3,sn(:,1:nlay))
634!        call writediagfi(ngrid,'vdif_kc_sm','','',3,sm(:,1:nlay))
635
[38]636      RETURN
637      END
Note: See TracBrowser for help on using the repository browser.