source: LMDZ4/trunk/libf/phytherm/yamada4.F @ 852

Last change on this file since 852 was 852, checked in by Laurent Fairhead, 17 years ago

Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007
LF

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