source: trunk/libf/phyvenus/yamada4.F @ 4

Last change on this file since 4 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

  • libf/phytitan : physique de Titan
  • libf/chimtitan: chimie de Titan
  • libf/phyvenus : physique de Venus
File size: 15.9 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/yamada4.F,v 1.1 2004/06/22 11:45:36 lmdzadmin Exp $
3!
4      SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
5c     s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
6     s   ,zlev,zlay,u,v,teta,cd,km,kn,kq,ustar
7     s   ,iflag_pbl)
8      IMPLICIT NONE
9c.......................................................................
10#include "dimensions.h"
11#include "dimphy.h"
12c.......................................................................
13c
14c dt : pas de temps
15c g  : g
16c zlev : altitude a chaque niveau (interface inferieure de la couche
17c        de meme indice)
18c zlay : altitude au centre de chaque couche
19c u,v : vitesse au centre de chaque couche
20c       (en entree : la valeur au debut du pas de temps)
21c teta : temperature potentielle au centre de chaque couche
22c        (en entree : la valeur au debut du pas de temps)
23c cd : cdrag
24c      (en entree : la valeur au debut du pas de temps)
25c q2 : $q^2$ au bas de chaque couche
26c      (en entree : la valeur au debut du pas de temps)
27c      (en sortie : la valeur a la fin du pas de temps)
28c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
29c      couche)
30c      (en sortie : la valeur a la fin du pas de temps)
31c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
32c      (en sortie : la valeur a la fin du pas de temps)
33c
34c  iflag_pbl doit valoir entre 6 et 9
35c      l=6, on prend  systematiquement une longueur d'equilibre
36c    iflag_pbl=6 : MY 2.0
37c    iflag_pbl=7 : MY 2.0.Fournier
38c    iflag_pbl=8 : MY 2.5
39c    iflag_pbl=9 : un test ?
40
41c.......................................................................
42      REAL dt,g,rconst
43      real plev(klon,klev+1),temp(klon,klev)
44      real ustar(klon)
45      real kmin,qmin,pblhmin(klon),coriol(klon)
46      REAL zlev(klon,klev+1)
47      REAL zlay(klon,klev)
48      REAL u(klon,klev)
49      REAL v(klon,klev)
50      REAL teta(klon,klev)
51      REAL cd(klon)
52      REAL q2(klon,klev+1),qpre
53      REAL unsdz(klon,klev)
54      REAL unsdzdec(klon,klev+1)
55
56      REAL km(klon,klev+1)
57      REAL kmpre(klon,klev+1),tmp2
58      REAL mpre(klon,klev+1)
59      REAL kn(klon,klev+1)
60      REAL kq(klon,klev+1)
61      real ff(klon,klev+1),delta(klon,klev+1)
62      real aa(klon,klev+1),aa0,aa1
63      integer iflag_pbl,ngrid
64
65
66      integer nlay,nlev
67      PARAMETER (nlay=klev)
68      PARAMETER (nlev=klev+1)
69
70      logical first
71      integer ipas
72      save first,ipas
73      data first,ipas/.true.,0/
74
75
76      integer ig,k
77
78
79      real ri,zrif,zalpha,zsm,zsn
80      real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
81
82      real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
83      real dtetadz(klon,klev+1)
84      real m2cstat,mcstat,kmcstat
85      real l(klon,klev+1),l0(klon)
86      save l0
87c  ATTENTION! mis ici car j'ai enlevé q2 des arguments...
88c   sinon, c'est au-dessus que ça se passe...
89      save q2
90
91      real sq(klon),sqz(klon),zz(klon,klev+1)
92      integer iter
93
94      real ric,rifc,b1,kap
95      save ric,rifc,b1,kap
96      data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
97
98      real frif,falpha,fsm
99      real fl,zzz,zl0,zq2,zn2
100
101      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
102     s  ,lyam(klon,klev),knyam(klon,klev)
103     s  ,w2yam(klon,klev),t2yam(klon,klev)
104      common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam
105
106      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
107      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
108      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
109      fl(zzz,zl0,zq2,zn2)=
110     s     max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
111     s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
112
113      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
114           stop'probleme de coherence dans appel a MY'
115      endif
116
117c===================================
118c INITIALISATIONS (surtout pour k=1, à cause diagnostiques...)
119      if (first) then
120        q2 = 0.
121      endif
122c===================================
123     
124      ipas=ipas+1
125      if (0.eq.1.and.first) then
126      do ig=1,1000
127         ri=(ig-800.)/500.
128         if (ri.lt.ric) then
129            zrif=frif(ri)
130         else
131            zrif=rifc
132         endif
133         if(zrif.lt.0.16) then
134            zalpha=falpha(zrif)
135            zsm=fsm(zrif)
136         else
137            zalpha=1.12
138            zsm=0.085
139         endif
140c     print*,ri,rif,zalpha,zsm
141      enddo
142      endif
143
144c.......................................................................
145c  les increments verticaux
146c.......................................................................
147c
148c!!!!! allerte !!!!!c
149c!!!!! zlev n'est pas declare a nlev !!!!!c
150c!!!!! ---->
151                                                      DO ig=1,ngrid
152            zlev(ig,nlev)=zlay(ig,nlay)
153     &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
154                                                      ENDDO
155c!!!!! <----
156c!!!!! allerte !!!!!c
157c
158      DO k=1,nlay
159                                                      DO ig=1,ngrid
160        unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
161                                                      ENDDO
162      ENDDO
163                                                      DO ig=1,ngrid
164      unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
165                                                      ENDDO
166      DO k=2,nlay
167                                                      DO ig=1,ngrid
168        unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
169                                                     ENDDO
170      ENDDO
171                                                      DO ig=1,ngrid
172      unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
173                                                     ENDDO
174c
175c.......................................................................
176
177c===================================
178c INITIALISATIONS (surtout pour k=1, à cause diagnostiques...)
179        dz = 0.
180        m2 = 0.
181        dtetadz = 0.
182        n2 = 0.
183        rif = 0.
184        alpha = 0.
185        sm = 0.
186        zz = 0.
187        l = 0.
188c===================================
189      do k=2,klev
190                                                          do ig=1,ngrid
191         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
192         m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
193     s             /(dz(ig,k)*dz(ig,k))
194         dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
195         n2(ig,k)=g*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
196c        n2(ig,k)=0.
197         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
198         if (ri.lt.ric) then
199            rif(ig,k)=frif(ri)
200         else
201            rif(ig,k)=rifc
202         endif
203         if(rif(ig,k).lt.0.16) then
204            alpha(ig,k)=falpha(rif(ig,k))
205            sm(ig,k)=fsm(rif(ig,k))
206         else
207            alpha(ig,k)=1.12
208            sm(ig,k)=0.085
209         endif
210         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
211c     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
212
213
214                                                          enddo
215      enddo
216
217
218c====================================================================
219c   Au premier appel, on determine l et q2 de facon iterative.
220c iterration pour determiner la longueur de melange
221
222
223      if (first.or.iflag_pbl.eq.6) then
224                                                          do ig=1,ngrid
225      l0(ig)=10.
226                                                          enddo
227      do k=2,klev-1
228                                                          do ig=1,ngrid
229        l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
230                                                          enddo
231      enddo
232
233      do iter=1,10
234                                                          do ig=1,ngrid
235         sq(ig)=1.e-10
236         sqz(ig)=1.e-10
237                                                          enddo
238         do k=2,klev-1
239                                                          do ig=1,ngrid
240           q2(ig,k)=l(ig,k)**2*zz(ig,k)
241           l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
242           zq=sqrt(q2(ig,k))
243           sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
244           sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
245                                                          enddo
246         enddo
247                                                          do ig=1,ngrid
248         l0(ig)=0.2*sqz(ig)/sq(ig)
249c        l0(ig)=30.
250                                                          enddo
251c      print*,'ITER=',iter,'  L0=',l0
252
253      enddo
254
255c     print*,'Fin de l initialisation de q2 et l0'
256
257      endif ! first
258
259c====================================================================
260c  Calcul de la longueur de melange.
261c====================================================================
262
263c   Mise a jour de l0
264                                                          do ig=1,ngrid
265      sq(ig)=1.e-10
266      sqz(ig)=1.e-10
267                                                          enddo
268      do k=2,klev-1
269                                                          do ig=1,ngrid
270        zq=sqrt(q2(ig,k))
271        sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
272        sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
273                                                          enddo
274      enddo
275                                                          do ig=1,ngrid
276      l0(ig)=0.2*sqz(ig)/sq(ig)
277c        l0(ig)=30.
278                                                          enddo
279c      print*,'ITER=',iter,'  L0=',l0
280c   calcul de l(z)
281      do k=2,klev
282                                                          do ig=1,ngrid
283         l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
284         if(first) then
285           q2(ig,k)=l(ig,k)**2*zz(ig,k)
286         endif
287                                                          enddo
288      enddo
289
290c====================================================================
291c   Yamada 2.0
292c====================================================================
293      if (iflag_pbl.eq.6) then
294
295      do k=2,klev
296                                                          do ig=1,ngrid
297         q2(ig,k)=l(ig,k)**2*zz(ig,k)
298                                                          enddo
299      enddo
300
301
302      else if (iflag_pbl.eq.7) then
303c====================================================================
304c   Yamada 2.Fournier
305c====================================================================
306
307c  Calcul de l,  km, au pas precedent
308      do k=2,klev
309                                                          do ig=1,ngrid
310c        print*,'SMML=',sm(ig,k),l(ig,k)
311         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
312         kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
313         mpre(ig,k)=sqrt(m2(ig,k))
314c        print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
315                                                          enddo
316      enddo
317
318      do k=2,klev-1
319                                                          do ig=1,ngrid
320        m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
321        mcstat=sqrt(m2cstat)
322
323c        print*,'M2 L=',k,mpre(ig,k),mcstat
324c
325c  -----{puis on ecrit la valeur de q qui annule l'equation de m
326c        supposee en q3}
327c
328        IF (k.eq.2) THEN
329          kmcstat=1.E+0 / mcstat
330     &    *( unsdz(ig,k)*kmpre(ig,k+1)
331     &                        *mpre(ig,k+1)
332     &      +unsdz(ig,k-1)
333     &              *cd(ig)
334     &              *( sqrt(u(ig,3)**2+v(ig,3)**2)
335     &                -mcstat/unsdzdec(ig,k)
336     &                -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
337     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
338        ELSE
339          kmcstat=1.E+0 / mcstat
340     &    *( unsdz(ig,k)*kmpre(ig,k+1)
341     &                        *mpre(ig,k+1)
342     &      +unsdz(ig,k-1)*kmpre(ig,k-1)
343     &                          *mpre(ig,k-1) )
344     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
345        ENDIF
346c       print*,'T2 L=',k,tmp2
347        tmp2=kmcstat
348     &      /( sm(ig,k)/q2(ig,k) )
349     &      /l(ig,k)
350        q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
351c       print*,'Q2 L=',k,q2(ig,k)
352c
353                                                          enddo
354      enddo
355
356      else if (iflag_pbl.ge.8) then
357c====================================================================
358c   Yamada 2.5 a la Didi
359c====================================================================
360
361
362c  Calcul de l,  km, au pas precedent
363      do k=2,klev
364                                                          do ig=1,ngrid
365c        print*,'SMML=',sm(ig,k),l(ig,k)
366         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
367         if (delta(ig,k).lt.1.e-20) then
368c     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
369            delta(ig,k)=1.e-20
370         endif
371         km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
372         aa0=
373     s   (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
374         aa1=
375     s   (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
376c abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
377         aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))
378c     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
379         qpre=sqrt(q2(ig,k))
380         if (iflag_pbl.eq.8 ) then
381            if (aa(ig,k).gt.0.) then
382               q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
383            else
384               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
385            endif
386         else ! iflag_pbl=9
387            if (aa(ig,k)*qpre.gt.0.9) then
388               q2(ig,k)=(qpre*10.)**2
389            else
390               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
391            endif
392         endif
393         q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
394c     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
395                                                          enddo
396      enddo
397
398      endif ! Fin du cas 8
399
400c     print*,'OK8'
401
402c====================================================================
403c   Calcul des coefficients de mélange
404c====================================================================
405      do k=2,klev
406c     print*,'k=',k
407                                                          do ig=1,ngrid
408cabde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
409         zq=sqrt(q2(ig,k))
410         km(ig,k)=l(ig,k)*zq*sm(ig,k)
411         kn(ig,k)=km(ig,k)*alpha(ig,k)
412         kq(ig,k)=l(ig,k)*zq*0.2
413c     print*,'KML=',km(ig,k),kn(ig,k)
414                                                          enddo
415      enddo
416
417c     if (iflag_pbl.ge.7..and.0.eq.1) then
418c        q2(:,1)=q2(:,2)
419c        call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
420c     endif
421
422c   Traitement des cas noctrunes avec l'introduction d'une longueur
423c   minilale.
424
425c====================================================================
426c   Traitement particulier pour les cas tres stables.
427c   D'apres Holtslag Boville.
428
429c     print*,'YAMADA4 0'
430
431                                                          do ig=1,ngrid
432      coriol(ig)=1.e-4
433      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
434                                                          enddo
435      if (first) then
436       print*,'A REVOIR!! coriol ?? pblhmin ',pblhmin
437      endif
438CTest a remettre 21 11 02
439c test abd 13 05 02      if(0.eq.1) then
440      if(1.eq.1) then
441      do k=2,klev
442         do ig=1,klon
443            if (teta(ig,2).gt.teta(ig,1)) then
444               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
445               kmin=kap*zlev(ig,k)*qmin
446            else
447               kmin=-1. ! kmin n'est utilise que pour les SL stables.
448            endif
449            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
450c               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
451c     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
452               kn(ig,k)=kmin
453               km(ig,k)=kmin
454               kq(ig,k)=kmin
455c   la longueur de melange est suposee etre l= kap z
456c   K=l q Sm d'ou q2=(K/l Sm)**2
457               q2(ig,k)=(qmin/sm(ig,k))**2
458            endif
459         enddo
460      enddo
461      endif
462
463c     print*,'YAMADA4 1'
464c   Diagnostique pour stockage
465
466      rino=rif
467      smyam(:,1:klev)=sm(:,1:klev)
468      styam=sm(:,1:klev)*alpha(:,1:klev)
469      lyam(1:klon,1:klev)=l(:,1:klev)
470      knyam(1:klon,1:klev)=kn(:,1:klev)
471
472c   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
473
474        if(1.eq.0)then
475      w2yam=q2(:,1:klev)*0.24
476     s    +lyam(:,1:klev)*5.17*kn(:,1:klev)*n2(:,1:klev)
477     s   /sqrt(q2(:,1:klev))
478
479      t2yam=9.1*kn(:,1:klev)*dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))
480     s  *lyam(:,1:klev)
481        endif
482
483c     print*,'OKFIN'
484      first=.false.
485      return
486      end
Note: See TracBrowser for help on using the repository browser.