source: trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F @ 3523

Last change on this file since 3523 was 3236, checked in by gmilcareck, 9 months ago

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