My Project
 All Classes Files Functions Variables Macros
yamada_c.F90
Go to the documentation of this file.
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 
120 REAL, DIMENSION(klon,klev+1) :: fluxu,fluxv,fluxt
121 REAL, DIMENSION(klon,klev+1) :: dddu,dddv,dddt
122 REAL, DIMENSION(klon,klev) :: exner,masse
123 REAL, 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
144 if (okiophys==1) then
145 call iophys_ecrit('q2i',klev,'q2 debut my','m2/s2',q2(:,1:klev))
146 call iophys_ecrit('kmi',klev,'Kz debut my','m/s2',km(:,1:klev))
147 endif
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
208 if (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))
222 endif
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
344 if (okiophys==1) then
345 call iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev))
346 call iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev))
347 call iophys_ecrit('Km2',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))
348 call iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev))
349 endif
350 #endif
351 
352 
353 IF (iflag_pbl<20) then
354  ! For diagnostics only
355  RETURN
356 
357 ELSE
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
439 if (okiophys==1) then
440 call iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev))
441 call iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev))
442 call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
443 call iophys_ecrit('q2neg',klev,'KTE non bornee','m2/s',q2neg(:,1:klev))
444 call iophys_ecrit('alpha',klev,'alpha','',alpha(:,1:klev))
445 call iophys_ecrit('sm',klev,'sm','',sm(:,1:klev))
446 call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
447 call iophys_ecrit('kmf',klev,'Kz final','m2/s',km(:,1:klev))
448 call iophys_ecrit('knf',klev,'Kz final','m2/s',kn(:,1:klev))
449 call iophys_ecrit('kqf',klev,'Kz final','m2/s',kq(:,1:klev))
450 endif
451 #endif
452 
453 ENDIF
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