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

Last change on this file since 2156 was 1779, checked in by aslmd, 7 years ago

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