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

Last change on this file since 3730 was 3662, checked in by gmilcareck, 4 months ago

Generic PCM:
Minor changing for the virtual potential temperature correction.
And convadj.F becomes convadj.F90 .
GM
File size: 19.6 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,
6     &                        virtual_theta_correction
7      IMPLICIT NONE
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      INTEGER,INTENT(IN) :: ngrid
32      INTEGER,INTENT(IN) :: nlay
33      INTEGER,INTENT(IN) :: nq
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) :: zq(ngrid,nlay,nq)
40      REAL,INTENT(IN) :: teta(ngrid,nlay)
41      REAL,INTENT(IN) :: cd(ngrid)
42      REAL,INTENT(INOUT) :: q2(ngrid,nlay+1)
43      REAL,INTENT(OUT) :: km(ngrid,nlay+1)
44      REAL,INTENT(OUT) :: kn(ngrid,nlay+1)
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 :: nlev
58      REAL unsdz(ngrid,nlay)
59      REAL unsdzdec(ngrid,nlay+1)
60      REAL q(ngrid,nlay+1)
61      REAL tetav(ngrid,nlay)
62      integer iq,igcm_generic_vap, igcm_generic_ice
63      logical call_ice_vap_generic
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,PARAMETER :: gnmin=-10.E+0
125      REAL,PARAMETER :: gnmax=0.0233E+0
126      LOGICAL gninf
127      LOGICAL gnsup
128      REAL gm
129c      REAL ri(ngrid,nlaye+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,PARAMETER :: kappa=0.4E+0
145      REAL,PARAMETER :: long0=160.E+0
146      REAL,PARAMETER :: a1=0.92E+0,a2=0.74E+0
147      REAL,PARAMETER :: b1=16.6E+0,b2=10.1E+0,c1=0.08E+0
148      REAL,PARAMETER :: cn1=a2*(1.E+0 -6.E+0 *a1/b1)
149      REAL,PARAMETER :: cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
150      REAL,PARAMETER :: cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
151      REAL,PARAMETER :: cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*
152     &          (1.E+0 -6.E+0 *a1/b1)-3.E+0 *c1*(b2+6.E+0 *a1)))
153      REAL,PARAMETER :: cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
154      REAL,PARAMETER :: cm4=-9.E+0 *a1*a2
155c.......................................................................
156c
157c termq : termes en $q$ dans l'equation de q2
158c termq3 : termes en $q^3$ dans l'equation de q2
159c termqm2 : termes en $q*m^2$ dans l'equation de q2
160c termq3m2 : termes en $q^3*m^2$ dans l'equation de q2
161c
162c.......................................................................
163      REAL termq
164      REAL termq3
165      REAL termqm2
166      REAL termq3m2
167c.......................................................................
168c
169c q2min : borne inferieure de q2
170c q2max : borne superieure de q2
171c
172c.......................................................................
173      REAL,PARAMETER :: q2min=1.E-3
174      REAL,PARAMETER :: q2max=1.E+2
175c.......................................................................
176c knmin : borne inferieure de kn
177c kmmin : borne inferieure de km
178c.......................................................................
179      REAL,PARAMETER :: knmin=1.E-5
180      REAL,PARAMETER :: kmmin=1.E-5
181c.......................................................................
182      INTEGER ilay,ilev,igrid
183      REAL tmp1,tmp2
184c.......................................................................
185c
186
187! initialization of local variables:
188      nlev=nlay+1
189      long(:,:)=0.
190      n2(:,:)=0.
191      sn(:,:)=0.
192      snq2(:,:)=0.
193      sm(:,:)=0.
194      smq2(:,:)=0.
195
196c.......................................................................
197c  traitment des valeur de q2 en entree
198c.......................................................................
199c
200      DO ilev=1,nlev
201       DO igrid=1,ngrid
202        q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
203        q(igrid,ilev)=sqrt(q2(igrid,ilev))
204       ENDDO
205      ENDDO
206c
207      DO igrid=1,ngrid
208        tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
209        q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
210        q2(igrid,1)=amax1(q2(igrid,1),q2min)
211        q(igrid,1)=sqrt(q2(igrid,1))
212      ENDDO
213c
214c.......................................................................
215c  les increments verticaux
216c.......................................................................
217c
218c!!!!! allerte !!!!!c
219c!!!!! zlev n'est pas declare a nlev !!!!!c
220c!!!!! ---->
221c                                                     DO igrid=1,ngrid
222c           zlev(igrid,nlev)=zlay(igrid,nlay)
223c    &             +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
224c                                                     ENDDO           
225c!!!!! <----
226c!!!!! allerte !!!!!c
227c
228      DO ilay=1,nlay
229       DO igrid=1,ngrid
230        unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
231       ENDDO
232      ENDDO
233
234      DO igrid=1,ngrid
235        unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
236      ENDDO
237
238      DO ilay=2,nlay
239        DO igrid=1,ngrid
240          unsdzdec(igrid,ilay)=1.E+0/
241     &                          (zlay(igrid,ilay)-zlay(igrid,ilay-1))
242        ENDDO
243      ENDDO
244     
245      DO igrid=1,ngrid
246        unsdzdec(igrid,nlay+1)=1.E+0/
247     &                          (zlev(igrid,nlay+1)-zlay(igrid,nlay))
248      ENDDO
249c
250c.......................................................................
251c  le cisaillement et le gradient de temperature
252c.......................................................................
253c
254      DO igrid=1,ngrid
255        m2(igrid,1)=(unsdzdec(igrid,1)
256     &                   *u(igrid,1))**2
257     &                 +(unsdzdec(igrid,1)
258     &                   *v(igrid,1))**2
259        m(igrid,1)=sqrt(m2(igrid,1))
260        mpre(igrid,1)=m(igrid,1)
261      ENDDO
262
263c
264c.......................................................................
265c  Virtual theta correction
266c.......................................................................
267c
268
269      if((generic_condensation) .and. (virtual_theta_correction)) THEN
270        DO ilev=1,nlay
271         DO igrid=1,ngrid
272          DO iq=1,nq
273            call generic_tracer_index(nq,iq,igcm_generic_vap,
274     &         igcm_generic_ice,call_ice_vap_generic)
275            if(call_ice_vap_generic) then
276               tetav(igrid,ilev) = teta(igrid,ilev)*
277     &               (1.e0+(1.e0/epsi_generic-1.e0)*
278     &               zq(igrid,ilev,igcm_generic_vap))
279            endif
280          ENDDO
281         ENDDO
282        ENDDO
283      endif
284
285c
286c-----------------------------------------------------------------------
287      DO ilev=2,nlev-1
288       DO igrid=1,ngrid
289c-----------------------------------------------------------------------
290c
291      if((generic_condensation) .and. (virtual_theta_correction)) THEN
292        n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
293     &                 *(tetav(igrid,ilev)-tetav(igrid,ilev-1))
294     &                 /(tetav(igrid,ilev)+tetav(igrid,ilev-1))*2.E+0
295      else
296        n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
297     &                 *(teta(igrid,ilev)-teta(igrid,ilev-1))
298     &                 /(teta(igrid,ilev)+teta(igrid,ilev-1)) *2.E+0
299      endif
300c
301c --->
302c       on ne sais traiter que les cas stratifies. et l'ajustement
303c       convectif est cense faire en sorte que seul des configurations
304c       stratifiees soient rencontrees en entree de cette routine.
305c       mais, bon ... on sait jamais (meme on sait que n2 prends
306c       quelques valeurs negatives ... parfois) alors :
307c<---
308c
309        IF (n2(igrid,ilev).lt.0.E+0) THEN
310          n2(igrid,ilev)=0.E+0
311        ENDIF
312c
313        m2(igrid,ilev)=(unsdzdec(igrid,ilev)
314     &                     *(u(igrid,ilev)-u(igrid,ilev-1)))**2
315     &                   +(unsdzdec(igrid,ilev)
316     &                     *(v(igrid,ilev)-v(igrid,ilev-1)))**2
317        m(igrid,ilev)=sqrt(m2(igrid,ilev))
318        mpre(igrid,ilev)=m(igrid,ilev)
319c
320c-----------------------------------------------------------------------
321       ENDDO
322      ENDDO
323c-----------------------------------------------------------------------
324c
325      DO igrid=1,ngrid
326        m2(igrid,nlev)=m2(igrid,nlev-1)
327        m(igrid,nlev)=m(igrid,nlev-1)
328        mpre(igrid,nlev)=m(igrid,nlev)
329      ENDDO
330     
331c
332c.......................................................................
333c  calcul des fonctions de stabilite
334c.......................................................................
335c
336c-----------------------------------------------------------------------
337      DO ilev=2,nlev-1
338       DO igrid=1,ngrid
339c-----------------------------------------------------------------------
340c
341        tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
342        long(igrid,ilev)=tmp1/(1.E+0 + tmp1/long0)
343        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
344     &                                           * n2(igrid,ilev)
345        gm=long(igrid,ilev)**2 / q2(igrid,ilev)
346     &                                           * m2(igrid,ilev)
347c
348        gninf=.false.
349        gnsup=.false.
350        long(igrid,ilev)=long(igrid,ilev)
351        long(igrid,ilev)=long(igrid,ilev)
352c
353        IF (gn.lt.gnmin) THEN
354          gninf=.true.
355          gn=gnmin
356        ENDIF
357c
358        IF (gn.gt.gnmax) THEN
359          gnsup=.true.
360          gn=gnmax
361        ENDIF
362c
363        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
364        sm(igrid,ilev)=
365     &    (cm1+cm2*gn)
366     &   /( (1.E+0 +cm3*gn)
367     &     *(1.E+0 +cm4*gn) )
368c
369        IF ((gninf).or.(gnsup)) THEN
370          snq2(igrid,ilev)=0.E+0
371          smq2(igrid,ilev)=0.E+0
372        ELSE
373          snq2(igrid,ilev)=
374     &     -gn
375     &     *(-cn1*cn2/(1.E+0 +cn2*gn)**2 )
376          smq2(igrid,ilev)=
377     &     -gn
378     &     *( cm2*(1.E+0 +cm3*gn)
379     &           *(1.E+0 +cm4*gn)
380     &       -( cm3*(1.E+0 +cm4*gn)
381     &         +cm4*(1.E+0 +cm3*gn) )
382     &       *(cm1+cm2*gn)            )
383     &     /( (1.E+0 +cm3*gn)
384     &       *(1.E+0 +cm4*gn) )**2
385        ENDIF
386c
387c --->
388c       la decomposition de Taylor en q2 n'a de sens que
389c       dans les cas stratifies ou sn et sm sont quasi
390c       proportionnels a q2. ailleurs on laisse le meme
391c       algorithme car l'ajustement convectif fait le travail.
392c       mais c'est delirant quand sn et snq2 n'ont pas le meme
393c       signe : dans ces cas, on ne fait pas la decomposition.
394c<---
395c
396        IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.E+0)
397     &      snq2(igrid,ilev)=0.E+0
398        IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.E+0)
399     &      smq2(igrid,ilev)=0.E+0
400c
401c-----------------------------------------------------------------------
402       ENDDO ! of DO igrid=1,ngrid
403      ENDDO ! of DO ilev=2,nlev-1
404c-----------------------------------------------------------------------
405c
406c.......................................................................
407c  calcul de km et kn au debut du pas de temps
408c.......................................................................
409c
410      DO igrid=1,ngrid
411        kn(igrid,1)=knmin
412        km(igrid,1)=kmmin
413        kmpre(igrid,1)=km(igrid,1)
414      ENDDO
415c
416c-----------------------------------------------------------------------
417      DO ilev=2,nlev-1
418       DO igrid=1,ngrid
419        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
420     &                                         *sn(igrid,ilev)
421        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
422     &                                         *sm(igrid,ilev)
423        kmpre(igrid,ilev)=km(igrid,ilev)
424       ENDDO
425      ENDDO
426c-----------------------------------------------------------------------
427c
428      DO igrid=1,ngrid
429        kn(igrid,nlev)=kn(igrid,nlev-1)
430        km(igrid,nlev)=km(igrid,nlev-1)
431        kmpre(igrid,nlev)=km(igrid,nlev)
432      ENDDO
433c
434c.......................................................................
435c  boucle sur les niveaux 2 a nlev-1
436c.......................................................................
437c
438c---->
439      DO 10001 ilev=2,nlev-1
440c---->
441      DO 10002 igrid=1,ngrid
442c
443c.......................................................................
444c
445c  calcul des termes sources et puits de l'equation de q2
446c  ------------------------------------------------------
447c
448        knq3=kn(igrid,ilev)*snq2(igrid,ilev)
449     &                                    /sn(igrid,ilev)
450        kmq3=km(igrid,ilev)*smq2(igrid,ilev)
451     &                                    /sm(igrid,ilev)
452c
453        termq=0.E+0
454        termq3=0.E+0
455        termqm2=0.E+0
456        termq3m2=0.E+0
457c
458        tmp1=dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
459        tmp2=dt*2.E+0 *kmq3*m2(igrid,ilev)
460        termqm2=termqm2
461     &    +dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
462     &    -dt*2.E+0 *kmq3*m2(igrid,ilev)
463        termq3m2=termq3m2
464     &    +dt*2.E+0 *kmq3*m2(igrid,ilev)
465c
466        termq=termq
467     &    -dt*2.E+0 *kn(igrid,ilev)*n2(igrid,ilev)
468     &    +dt*2.E+0 *knq3*n2(igrid,ilev)
469        termq3=termq3
470     &    -dt*2.E+0 *knq3*n2(igrid,ilev)
471c
472        termq3=termq3
473     &    -dt*2.E+0 *q(igrid,ilev)**3 / (b1*long(igrid,ilev))
474c
475c.......................................................................
476c
477c  resolution stationnaire couplee avec le gradient de vitesse local
478c  -----------------------------------------------------------------
479c
480c  -----{on cherche le cisaillement qui annule l'equation de q^2
481c        supposee en q3}
482c
483        tmp1=termq+termq3
484        tmp2=termqm2+termq3m2
485        m2cstat=m2(igrid,ilev)
486     &      -(tmp1+tmp2)/(dt*2.E+0*km(igrid,ilev))
487        mcstat=sqrt(m2cstat)
488c
489c  -----{puis on ecrit la valeur de q qui annule l'equation de m
490c        supposee en q3}
491c
492        IF (ilev.eq.2) THEN
493          kmcstat=1.E+0 / mcstat
494     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
495     &                        *mpre(igrid,ilev+1)
496     &      +unsdz(igrid,ilev-1)
497     &              *cd(igrid)
498     &              *( sqrt(u(igrid,3)**2+v(igrid,3)**2)
499     &                -mcstat/unsdzdec(igrid,ilev)
500     &                -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
501     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
502        ELSE
503          kmcstat=1.E+0 / mcstat
504     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
505     &                        *mpre(igrid,ilev+1)
506     &      +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
507     &                          *mpre(igrid,ilev-1) )
508     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
509        ENDIF
510        tmp2=kmcstat
511     &      /( sm(igrid,ilev)/q2(igrid,ilev) )
512     &      /long(igrid,ilev)
513        qcstat=tmp2**(1.E+0/3.E+0)
514        q2cstat=qcstat**2
515c
516c.......................................................................
517c
518c  choix de la solution finale
519c  ---------------------------
520c
521          q(igrid,ilev)=qcstat
522          q2(igrid,ilev)=q2cstat
523          m(igrid,ilev)=mcstat
524          m2(igrid,ilev)=m2cstat
525c
526c --->
527c       pour des raisons simples q2 est minore
528c<---
529c
530        IF (q2(igrid,ilev).lt.q2min) THEN
531          q2(igrid,ilev)=q2min
532          q(igrid,ilev)=sqrt(q2min)
533        ENDIF
534c
535c.......................................................................
536c
537c  calcul final de kn et km
538c  ------------------------
539c
540        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
541     &                                           * n2(igrid,ilev)
542        IF (gn.lt.gnmin) gn=gnmin
543        IF (gn.gt.gnmax) gn=gnmax
544        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
545        sm(igrid,ilev)=
546     &    (cm1+cm2*gn)
547     &   /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
548        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
549     &                 *sn(igrid,ilev)
550        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
551     &                 *sm(igrid,ilev)
552c
553c.......................................................................
554c
55510002 CONTINUE
556c
55710001 CONTINUE
558c
559c.......................................................................
560c
561c
562      DO igrid=1,ngrid
563        kn(igrid,1)=knmin
564        km(igrid,1)=kmmin
565        q2(igrid,nlev)=q2(igrid,nlev-1)
566        q(igrid,nlev)=q(igrid,nlev-1)
567        kn(igrid,nlev)=kn(igrid,nlev-1)
568        km(igrid,nlev)=km(igrid,nlev-1)
569      ENDDO
570
571      END
Note: See TracBrowser for help on using the repository browser.