source: trunk/LMDZ.PLUTO/libf/phypluto/vdif_kc.F @ 3380

Last change on this file since 3380 was 3184, checked in by afalco, 10 months ago

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