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

Last change on this file since 3948 was 3727, checked in by emillour, 7 months ago

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