source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/yamada4.F @ 3400

Last change on this file since 3400 was 634, checked in by Laurent Fairhead, 20 years ago

Modifications faites à la physique pour la rendre parallele YM
Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester
les modifs pleinement avant leurs inclusions dans le tronc principal
LF

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