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