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

Last change on this file since 1292 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

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