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