source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/vdif_kc.f @ 224

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