source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/yamada_c.F90 @ 3069

Last change on this file since 3069 was 1910, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1860:1909 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 22.9 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE yamada_c(ngrid,timestep,plev,play &
5     &   ,pu,pv,pt,d_u,d_v,d_t,cd,q2,km,kn,kq,d_t_diss,ustar &
6     &   ,iflag_pbl,okiophys)
7      use dimphy
8      IMPLICIT NONE
9#include "iniprint.h"
10#include "YOMCST.h"
11!.......................................................................
12!ym#include "dimensions.h"
13!ym#include "dimphy.h"
14!.......................................................................
15!
16! timestep : pas de temps
17! g  : g
18! zlev : altitude a chaque niveau (interface inferieure de la couche
19!        de meme indice)
20! zlay : altitude au centre de chaque couche
21! u,v : vitesse au centre de chaque couche
22!       (en entree : la valeur au debut du pas de temps)
23! teta : temperature potentielle au centre de chaque couche
24!        (en entree : la valeur au debut du pas de temps)
25! cd : cdrag
26!      (en entree : la valeur au debut du pas de temps)
27! q2 : $q^2$ au bas de chaque couche
28!      (en entree : la valeur au debut du pas de temps)
29!      (en sortie : la valeur a la fin du pas de temps)
30! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
31!      couche)
32!      (en sortie : la valeur a la fin du pas de temps)
33! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
34!      (en sortie : la valeur a la fin du pas de temps)
35!
36!  iflag_pbl doit valoir entre 6 et 9
37!      l=6, on prend  systematiquement une longueur d'equilibre
38!    iflag_pbl=6 : MY 2.0
39!    iflag_pbl=7 : MY 2.0.Fournier
40!    iflag_pbl=8/9 : MY 2.5
41!       iflag_pbl=8 with special obsolete treatments for convergence
42!       with Cmpi5 NPv3.1 simulations
43!    iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
44!    iflag_pbl=12 = 11 with vertical diffusion off q2
45!
46!  2013/04/01 (FH hourdin@lmd.jussieu.fr)
47!     Correction for very stable PBLs (iflag_pbl=10 and 11)
48!     iflag_pbl=8 converges numerically with NPv3.1
49!     iflag_pbl=11 -> the model starts with NP from start files created by ce0l
50!                  -> the model can run with longer time-steps.
51!.......................................................................
52
53      REAL, DIMENSION(klon,klev) :: d_u,d_v,d_t
54      REAL, DIMENSION(klon,klev) :: pu,pv,pt
55      REAL, DIMENSION(klon,klev) :: d_t_diss
56      INTEGER okiophys
57
58      REAL timestep
59      real plev(klon,klev+1)
60      real play(klon,klev)
61      real ustar(klon)
62      real kmin,qmin,pblhmin(klon),coriol(klon)
63      REAL zlev(klon,klev+1)
64      REAL zlay(klon,klev)
65      REAL zu(klon,klev)
66      REAL zv(klon,klev)
67      REAL zt(klon,klev)
68      REAL teta(klon,klev)
69      REAL cd(klon)
70      REAL q2(klon,klev+1),qpre
71      REAL unsdz(klon,klev)
72      REAL unsdzdec(klon,klev+1)
73
74      REAL km(klon,klev+1)
75      REAL kmpre(klon,klev+1),tmp2
76      REAL mpre(klon,klev+1)
77      REAL kn(klon,klev+1)
78      REAL kq(klon,klev+1)
79      real ff(klon,klev+1),delta(klon,klev+1)
80      real aa(klon,klev+1),aa0,aa1
81      integer iflag_pbl,ngrid
82      integer nlay,nlev
83
84      logical first
85      integer ipas
86      save first,ipas
87!FH/IM     data first,ipas/.true.,0/
88      data first,ipas/.false.,0/
89!$OMP THREADPRIVATE( first,ipas)
90
91      integer ig,k
92
93
94      real ri,zrif,zalpha,zsm,zsn
95      real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
96
97      real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
98      REAL, DIMENSION(klon,klev+1) :: km2,kn2,sqrtq
99      real dtetadz(klon,klev+1)
100      real m2cstat,mcstat,kmcstat
101      real l(klon,klev+1)
102      real leff(klon,klev+1)
103      real,allocatable,save :: l0(:)
104!$OMP THREADPRIVATE(l0)     
105      real sq(klon),sqz(klon),zz(klon,klev+1)
106      integer iter
107
108      real ric,rifc,b1,kap
109      save ric,rifc,b1,kap
110      data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
111!$OMP THREADPRIVATE(ric,rifc,b1,kap)
112      real frif,falpha,fsm
113      real fl,zzz,zl0,zq2,zn2
114
115      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
116      real lyam(klon,klev),knyam(klon,klev)
117      real w2yam(klon,klev),t2yam(klon,klev)
118      logical,save :: firstcall=.true.
119
120REAL, DIMENSION(klon,klev+1) :: fluxu,fluxv,fluxt
121REAL, DIMENSION(klon,klev+1) :: dddu,dddv,dddt
122REAL, DIMENSION(klon,klev) :: exner,masse
123REAL, DIMENSION(klon,klev+1) :: masseb,q2old,q2neg
124
125!$OMP THREADPRIVATE(firstcall)       
126      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
127      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
128      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
129      fl(zzz,zl0,zq2,zn2)= &
130     &     max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig)) &
131     &     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
132
133
134      if (firstcall) then
135        allocate(l0(klon))
136#ifdef IOPHYS
137        call iophys_ini
138#endif
139        firstcall=.false.
140      endif
141
142
143#ifdef IOPHYS
144if (okiophys==1) then
145call iophys_ecrit('q2i',klev,'q2 debut my','m2/s2',q2(:,1:klev))
146call iophys_ecrit('kmi',klev,'Kz debut my','m/s2',km(:,1:klev))
147endif
148#endif
149
150      nlay=klev
151      nlev=klev+1
152
153!-------------------------------------------------------------------------
154! Computation of conservative source terms from the turbulent tendencies
155!-------------------------------------------------------------------------
156
157
158   zu(:,:)=pu(:,:)+0.5*d_u(:,:)
159   zv(:,:)=pv(:,:)+0.5*d_v(:,:)
160   zt(:,:)=pt(:,:)+0.5*d_t(:,:)
161   do k=1,klev
162      exner(:,k)=(play(:,k)/plev(:,1))**RKAPPA
163      masse(:,k)=(plev(:,k)-plev(:,k+1))/RG
164   enddo
165   teta(:,:)=zt(:,:)/exner(:,:)
166
167! Atmospheric mass at layer interfaces, where the TKE is computed
168   masseb(:,:)=0.
169   do k=1,klev
170      masseb(:,k)=masseb(:,k)+masse(:,k)
171      masseb(:,k+1)=masseb(:,k+1)+masse(:,k)
172    enddo
173    masseb(:,:)=0.5*masseb(:,:)
174
175
176
177   zlev(:,1)=0.
178   zlay(:,1)=RCPD*teta(:,1)*(1.-exner(:,1))
179   do k=1,klev-1
180      zlay(:,k+1)=zlay(:,k)+0.5*RCPD*(teta(:,k)+teta(:,k+1))*(exner(:,k)-exner(:,k+1))/RG
181      zlev(:,k)=0.5*(zlay(:,k)+zlay(:,k+1)) ! PASBO
182   enddo
183
184   fluxu(:,klev+1)=0.
185   fluxv(:,klev+1)=0.
186   fluxt(:,klev+1)=0.
187
188   do k=klev,1,-1
189      fluxu(:,k)=fluxu(:,k+1)+masse(:,k)*d_u(:,k)
190      fluxv(:,k)=fluxv(:,k+1)+masse(:,k)*d_v(:,k)
191      fluxt(:,k)=fluxt(:,k+1)+masse(:,k)*d_t(:,k)/exner(:,k) ! Flux de theta
192   enddo
193
194   dddu(:,1)=2*zu(:,1)*fluxu(:,1)
195   dddv(:,1)=2*zv(:,1)*fluxv(:,1)
196   dddt(:,1)=(exner(:,1)-1.)*fluxt(:,1)
197
198   do k=2,klev
199      dddu(:,k)=(zu(:,k)-zu(:,k-1))*fluxu(:,k)
200      dddv(:,k)=(zv(:,k)-zv(:,k-1))*fluxv(:,k)
201      dddt(:,k)=(exner(:,k)-exner(:,k-1))*fluxt(:,k)
202   enddo
203   dddu(:,klev+1)=0.
204   dddv(:,klev+1)=0.
205   dddt(:,klev+1)=0.
206
207#ifdef IOPHYS
208if (okiophys==1) then
209      call iophys_ecrit('zlay',klev,'Geop','m',zlay)
210      call iophys_ecrit('teta',klev,'teta','K',teta)
211      call iophys_ecrit('temp',klev,'temp','K',zt)
212      call iophys_ecrit('pt',klev,'temp','K',pt)
213      call iophys_ecrit('d_u',klev,'d_u','m/s2',d_u)
214      call iophys_ecrit('d_v',klev,'d_v','m/s2',d_v)
215      call iophys_ecrit('d_t',klev,'d_t','K/s',d_t)
216      call iophys_ecrit('exner',klev,'exner','',exner)
217      call iophys_ecrit('masse',klev,'masse','',masse)
218      call iophys_ecrit('masseb',klev,'masseb','',masseb)
219      call iophys_ecrit('Cm2',klev,'m2 conserv','m/s',(dddu(:,1:klev)+dddv(:,1:klev))/(masseb(:,1:klev)*timestep))
220      call iophys_ecrit('Cn2',klev,'m2 conserv','m/s',(rcpd*dddt(:,1:klev)/masseb(:,1:klev))/timestep)
221      call iophys_ecrit('rifc',klev,'rif conservative','',rcpd*dddt(:,1:klev)/min(dddu(:,1:klev)+dddv(:,1:klev),-1.e-20))
222endif
223#endif
224
225
226
227      ipas=ipas+1
228
229
230!.......................................................................
231!  les increments verticaux
232!.......................................................................
233!
234!!!!!! allerte !!!!!c
235!!!!!! zlev n'est pas declare a nlev !!!!!c
236!!!!!! ---->
237                                                      DO ig=1,ngrid
238            zlev(ig,nlev)=zlay(ig,nlay) &
239     &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
240                                                      ENDDO
241!!!!!! <----
242!!!!!! allerte !!!!!c
243!
244      DO k=1,nlay
245                                                      DO ig=1,ngrid
246        unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
247                                                      ENDDO
248      ENDDO
249                                                      DO ig=1,ngrid
250      unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
251                                                      ENDDO
252      DO k=2,nlay
253                                                      DO ig=1,ngrid
254        unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
255                                                     ENDDO
256      ENDDO
257                                                      DO ig=1,ngrid
258      unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
259                                                     ENDDO
260!
261!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
262! Computing M^2, N^2, Richardson numbers, stability functions
263!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
264
265
266      do k=2,klev
267                                                          do ig=1,ngrid
268         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
269         m2(ig,k)=((zu(ig,k)-zu(ig,k-1))**2+(zv(ig,k)-zv(ig,k-1))**2)/(dz(ig,k)*dz(ig,k))
270         dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
271         n2(ig,k)=RG*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
272!        n2(ig,k)=0.
273         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
274         if (ri.lt.ric) then
275            rif(ig,k)=frif(ri)
276         else
277            rif(ig,k)=rifc
278         endif
279         if(rif(ig,k).lt.0.16) then
280            alpha(ig,k)=falpha(rif(ig,k))
281            sm(ig,k)=fsm(rif(ig,k))
282         else
283            alpha(ig,k)=1.12
284            sm(ig,k)=0.085
285         endif
286         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
287                                                          enddo
288      enddo
289
290
291
292!====================================================================
293!  Computing the mixing length
294!====================================================================
295
296!   Mise a jour de l0
297      if (iflag_pbl==8.or.iflag_pbl==10) then
298
299!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
300! Iterative computation of l0
301! This version is kept for iflag_pbl only for convergence
302! with NPv3.1 Cmip5 simulations
303!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
304
305                                                          do ig=1,ngrid
306      sq(ig)=1.e-10
307      sqz(ig)=1.e-10
308                                                          enddo
309      do k=2,klev-1
310                                                          do ig=1,ngrid
311        zq=sqrt(q2(ig,k))
312        sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
313        sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
314                                                          enddo
315      enddo
316                                                          do ig=1,ngrid
317      l0(ig)=0.2*sqz(ig)/sq(ig)
318                                                          enddo
319      do k=2,klev
320                                                          do ig=1,ngrid
321         l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
322                                                          enddo
323      enddo
324!     print*,'L0 cas 8 ou 10 ',l0
325
326      else
327
328!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329! In all other case, the assymptotic mixing length l0 is imposed (100m)
330!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
331
332          l0(:)=150.
333          do k=2,klev
334                                                          do ig=1,ngrid
335             l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
336                                                          enddo
337          enddo
338!     print*,'L0 cas autres ',l0
339
340      endif
341
342
343#ifdef IOPHYS
344if (okiophys==1) then
345call iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev))
346call iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev))
347call iophys_ecrit('Km2',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))
348call iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev))
349endif
350#endif
351
352
353IF (iflag_pbl<20) then
354      ! For diagnostics only
355      RETURN
356
357ELSE
358
359!  print*,'OK1'
360
361! Evolution of TKE under source terms K M2 and K N2
362   leff(:,:)=max(l(:,:),1.)
363   IF (iflag_pbl==29) THEN
364      km2(:,:)=km(:,:)*m2(:,:)
365      kn2(:,:)=kn2(:,:)*rif(:,:)
366   ELSEIF (iflag_pbl==25) THEN
367      DO k=1,klev
368         km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) &
369         &        /(masse(:,k)*timestep)
370         kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep)
371         leff(:,k)=0.5*(leff(:,k)+leff(:,k+1))
372      ENDDO
373      km2(:,klev+1)=0. ; kn2(:,klev+1)=0.
374   ELSE
375      km2(:,:)=-(dddu(:,:)+dddv(:,:))/(masseb(:,:)*timestep)
376      kn2(:,:)=rcpd*dddt(:,:)/(masseb(:,:)*timestep)
377   ENDIF
378   q2neg(:,:)=q2(:,:)+timestep*(km2(:,:)-kn2(:,:))
379   q2(:,:)=min(max(q2neg(:,:),1.e-10),1.e4)
380
381! Dissipation of TKE
382   q2old(:,:)=q2(:,:)
383   q2(:,:)=1./(1./sqrt(q2(:,:))+timestep/(2*leff(:,:)*b1))
384   q2(:,:)=q2(:,:)*q2(:,:)
385   IF (iflag_pbl<=24) THEN
386      DO k=1,klev
387         d_t_diss(:,k)=(masseb(:,k)*(q2neg(:,k)-q2(:,k))+masseb(:,k+1)*(q2neg(:,k+1)-q2(:,k+1)))/(2.*rcpd*masse(:,k))
388      ENDDO
389   ELSE IF (iflag_pbl<=27) THEN
390      DO k=1,klev
391         d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd
392      ENDDO
393   ENDIF
394!  print*,'iflag_pbl ',d_t_diss
395
396
397! Compuation of stability functions
398   IF (iflag_pbl/=29) THEN
399      DO k=1,klev
400      DO ig=1,ngrid
401         IF (ABS(km2(ig,k))<=1.e-20) THEN
402            rif(ig,k)=0.
403         ELSE
404            rif(ig,k)=min(kn2(ig,k)/km2(ig,k),rifc)
405         ENDIF
406         IF (rif(ig,k).lt.0.16) THEN
407            alpha(ig,k)=falpha(rif(ig,k))
408            sm(ig,k)=fsm(rif(ig,k))
409         else
410            alpha(ig,k)=1.12
411            sm(ig,k)=0.085
412         endif
413      ENDDO
414      ENDDO
415    ENDIF
416
417! Computation of turbulent diffusivities
418   IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN
419     DO k=2,klev
420        sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1)))
421     ENDDO
422   ELSE
423     DO k=2,klev
424        sqrtq(:,k)=sqrt(q2(:,k))
425     ENDDO
426   ENDIF
427   DO k=2,klev
428   DO ig=1,ngrid
429      km(ig,k)=leff(ig,k)*sqrtq(ig,k)*sm(ig,k)
430      kn(ig,k)=km(ig,k)*alpha(ig,k)
431      kq(ig,k)=leff(ig,k)*zq*0.2
432!         print*,q2(ig,k),zq,km(ig,k)
433   ENDDO
434   ENDDO
435
436
437
438#ifdef IOPHYS
439if (okiophys==1) then
440call iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev))
441call iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev))
442call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
443call iophys_ecrit('q2neg',klev,'KTE non bornee','m2/s',q2neg(:,1:klev))
444call iophys_ecrit('alpha',klev,'alpha','',alpha(:,1:klev))
445call iophys_ecrit('sm',klev,'sm','',sm(:,1:klev))
446call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
447call iophys_ecrit('kmf',klev,'Kz final','m2/s',km(:,1:klev))
448call iophys_ecrit('knf',klev,'Kz final','m2/s',kn(:,1:klev))
449call iophys_ecrit('kqf',klev,'Kz final','m2/s',kq(:,1:klev))
450endif
451#endif
452
453ENDIF
454
455
456!  print*,'OK2'
457      RETURN
458!====================================================================
459!   Yamada 2.0
460!====================================================================
461      if (iflag_pbl.eq.6) then
462
463      do k=2,klev
464         q2(:,k)=l(:,k)**2*zz(:,k)
465      enddo
466
467
468      else if (iflag_pbl.eq.7) then
469!====================================================================
470!   Yamada 2.Fournier
471!====================================================================
472
473!  Calcul de l,  km, au pas precedent
474      do k=2,klev
475                                                          do ig=1,ngrid
476!        print*,'SMML=',sm(ig,k),l(ig,k)
477         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
478         kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
479         mpre(ig,k)=sqrt(m2(ig,k))
480!        print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
481                                                          enddo
482      enddo
483
484      do k=2,klev-1
485                                                          do ig=1,ngrid
486        m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
487        mcstat=sqrt(m2cstat)
488
489!        print*,'M2 L=',k,mpre(ig,k),mcstat
490!
491!  -----{puis on ecrit la valeur de q qui annule l'equation de m
492!        supposee en q3}
493!
494        IF (k.eq.2) THEN
495          kmcstat=1.E+0 / mcstat &
496     &    *( unsdz(ig,k)*kmpre(ig,k+1) &
497     &                        *mpre(ig,k+1) &
498     &      +unsdz(ig,k-1) &
499     &              *cd(ig) &
500     &              *( sqrt(zu(ig,3)**2+zv(ig,3)**2) &
501     &                -mcstat/unsdzdec(ig,k) &
502     &                -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2) &
503     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
504        ELSE
505          kmcstat=1.E+0 / mcstat &
506     &    *( unsdz(ig,k)*kmpre(ig,k+1) &
507     &                        *mpre(ig,k+1) &
508     &      +unsdz(ig,k-1)*kmpre(ig,k-1) &
509     &                          *mpre(ig,k-1) ) &
510     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
511        ENDIF
512!       print*,'T2 L=',k,tmp2
513        tmp2=kmcstat &
514     &      /( sm(ig,k)/q2(ig,k) ) &
515     &      /l(ig,k)
516        q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
517!       print*,'Q2 L=',k,q2(ig,k)
518!
519                                                          enddo
520      enddo
521
522      else if (iflag_pbl==8.or.iflag_pbl==9) then
523!====================================================================
524!   Yamada 2.5 a la Didi
525!====================================================================
526
527
528!  Calcul de l,  km, au pas precedent
529      do k=2,klev
530                                                          do ig=1,ngrid
531!        print*,'SMML=',sm(ig,k),l(ig,k)
532         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
533         if (delta(ig,k).lt.1.e-20) then
534!     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
535            delta(ig,k)=1.e-20
536         endif
537         km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
538         aa0=(m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
539         aa1=(m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
540! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
541         aa(ig,k)=aa1*timestep/(delta(ig,k)*l(ig,k))
542!     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
543         qpre=sqrt(q2(ig,k))
544!        if (iflag_pbl.eq.8 ) then
545            if (aa(ig,k).gt.0.) then
546               q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
547            else
548               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
549            endif
550!        else ! iflag_pbl=9
551!           if (aa(ig,k)*qpre.gt.0.9) then
552!              q2(ig,k)=(qpre*10.)**2
553!           else
554!              q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
555!           endif
556!        endif
557         q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
558!     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
559                                                          enddo
560      enddo
561
562      else if (iflag_pbl>=10) then
563
564!        print*,'Schema mixte D'
565!        print*,'Longueur ',l(:,:)
566         do k=2,klev-1
567            l(:,k)=max(l(:,k),1.)
568            km(:,k)=l(:,k)*sqrt(q2(:,k))*sm(:,k)
569            q2(:,k)=q2(:,k)+timestep*km(:,k)*m2(:,k)*(1.-rif(:,k))
570            q2(:,k)=min(max(q2(:,k),1.e-10),1.e4)
571            q2(:,k)=1./(1./sqrt(q2(:,k))+timestep/(2*l(:,k)*b1))
572            q2(:,k)=q2(:,k)*q2(:,k)
573         enddo
574
575
576      else
577         stop'Cas nom prevu dans yamada4'
578
579      endif ! Fin du cas 8
580
581!     print*,'OK8'
582
583!====================================================================
584!   Calcul des coefficients de m�ange
585!====================================================================
586      do k=2,klev
587!     print*,'k=',k
588                                                          do ig=1,ngrid
589!abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
590         zq=sqrt(q2(ig,k))
591         km(ig,k)=l(ig,k)*zq*sm(ig,k)
592         kn(ig,k)=km(ig,k)*alpha(ig,k)
593         kq(ig,k)=l(ig,k)*zq*0.2
594!     print*,'KML=',km(ig,k),kn(ig,k)
595                                                          enddo
596      enddo
597
598! Transport diffusif vertical de la TKE.
599      if (iflag_pbl.ge.12) then
600!       print*,'YAMADA VDIF'
601        q2(:,1)=q2(:,2)
602        call vdif_q2(timestep,RG,RD,ngrid,plev,zt,kq,q2)
603      endif
604
605!   Traitement des cas noctrunes avec l'introduction d'une longueur
606!   minilale.
607
608!====================================================================
609!   Traitement particulier pour les cas tres stables.
610!   D'apres Holtslag Boville.
611
612      if (prt_level>1) THEN
613       print*,'YAMADA4 0'
614      endif !(prt_level>1) THEN
615                                                          do ig=1,ngrid
616      coriol(ig)=1.e-4
617      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
618                                                          enddo
619
620!     print*,'pblhmin ',pblhmin
621!Test a remettre 21 11 02
622! test abd 13 05 02      if(0.eq.1) then
623      if(1==1) then
624      if(iflag_pbl==8.or.iflag_pbl==10) then
625
626      do k=2,klev
627         do ig=1,ngrid
628            if (teta(ig,2).gt.teta(ig,1)) then
629               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
630               kmin=kap*zlev(ig,k)*qmin
631            else
632               kmin=-1. ! kmin n'est utilise que pour les SL stables.
633            endif
634            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
635!               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
636!     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
637               kn(ig,k)=kmin
638               km(ig,k)=kmin
639               kq(ig,k)=kmin
640!   la longueur de melange est suposee etre l= kap z
641!   K=l q Sm d'ou q2=(K/l Sm)**2
642               q2(ig,k)=(qmin/sm(ig,k))**2
643            endif
644         enddo
645      enddo
646
647      else
648
649      do k=2,klev
650         do ig=1,ngrid
651            if (teta(ig,2).gt.teta(ig,1)) then
652               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
653               kmin=kap*zlev(ig,k)*qmin
654            else
655               kmin=-1. ! kmin n'est utilise que pour les SL stables.
656            endif
657            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
658!               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
659!     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
660               kn(ig,k)=kmin
661               km(ig,k)=kmin
662               kq(ig,k)=kmin
663!   la longueur de melange est suposee etre l= kap z
664!   K=l q Sm d'ou q2=(K/l Sm)**2
665               sm(ig,k)=1.
666               alpha(ig,k)=1.
667               q2(ig,k)=min((qmin/sm(ig,k))**2,10.)
668               zq=sqrt(q2(ig,k))
669               km(ig,k)=l(ig,k)*zq*sm(ig,k)
670               kn(ig,k)=km(ig,k)*alpha(ig,k)
671               kq(ig,k)=l(ig,k)*zq*0.2
672            endif
673         enddo
674      enddo
675      endif
676
677      endif
678
679      if (prt_level>1) THEN
680       print*,'YAMADA4 1'
681      endif !(prt_level>1) THEN
682!   Diagnostique pour stokage
683
684      if(1.eq.0)then
685      rino=rif
686      smyam(1:ngrid,1)=0.
687      styam(1:ngrid,1)=0.
688      lyam(1:ngrid,1)=0.
689      knyam(1:ngrid,1)=0.
690      w2yam(1:ngrid,1)=0.
691      t2yam(1:ngrid,1)=0.
692
693      smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)
694      styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev)
695      lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev)
696      knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev)
697
698!   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
699
700      w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24 &
701     &    +lyam(1:ngrid,2:klev)*5.17*kn(1:ngrid,2:klev) &
702     &    *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev))
703
704      t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev) &
705     &    *dtetadz(1:ngrid,2:klev)**2 &
706     &    /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev)
707      endif
708
709!     print*,'OKFIN'
710      first=.false.
711      return
712      end
Note: See TracBrowser for help on using the repository browser.