source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phylmd/yamada_c.F90 @ 5005

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

New version of Mellor et Yamada pronostic TKE

... based on energy transfer from the mean state.

The new version is yamada_c.
It must be called after vertical diffusion rather than just before
since the source terms u_z w'u' + v_z w'theta' and g/theta w'theta'
are diagnosed from the vertical diffusion (as energy loss from the mean
state) rather than computed as K (u_z2+v_z2) or g/\theta K theta_z
(where _z means vertical derivative).
The call to this version is controled by iflag_pbl.

iflag_pbl = 20 : with TKE computed at interfaces between layers
iflag_pbl = 25 : with TKE computed within the layer
In both cases, the dissipation of turbulence is translated into heat, and
passed to the physics as dtdiss (temperature tendency due to dissipation

of TKE).

The diffusion coefficient being computed after dissipation, it must be
kept for diffusion at the next time step, and thus be stored in the
restartphy file.
This coefficient must be computed and stored for each sub-surface.

A new way of managing subsurface variables is introduced.
For arrays of the form X(:,nbsrf) are extented to X(:,nbsrf+1), where
is_ave=nbsrf+1, is an additional sub-surface which contains the averaged values.

coef_diff_turb_mod.F90 : change of flags.
ener_conserv.F90 : energy conservation must not be applied in those

cases

indicesol.h : definition of is_ave
pbl_surface_mod.F90 : call to yamada_c and changes in the management of

coefh/coefm

physiq.F : Change in the initialisation of pmflxr/s and

modified calls to pbl_surface_mod (introduction
of dtdiss, initialisation
of pbl_tke and coefh in 1D).

phys_local_var_mod.F90 : declaration of d_t_diss
phys_output_mod.F90 : new outputs (bils_tke and bils_diss) and

coefh->coefh(:,:,is_ave)

phys_state_var_mod.F90 : modified declaration for coefh and coefm

(nbsrf -> nbsrf+1)

FH

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.