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

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

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

  • Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
  • dyn3d/control.h is now module control_mod.F90
  • rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
  • added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
  • created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

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