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

Last change on this file since 1351 was 1308, checked in by emillour, 10 years ago

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