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

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

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