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