source: LMDZ4/trunk/libf/phylmd/yamada4.F @ 979

Last change on this file since 979 was 938, checked in by lmdzadmin, 17 years ago

Enleve prints par defaut
IM

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