source: LMDZ5/trunk/libf/phylmd/yamada_c.F90 @ 2373

Last change on this file since 2373 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 22.8 KB
RevLine 
[1761]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
[2311]8      USE print_control_mod, ONLY: prt_level
[1761]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
116REAL, DIMENSION(klon,klev+1) :: fluxu,fluxv,fluxt
117REAL, DIMENSION(klon,klev+1) :: dddu,dddv,dddt
118REAL, DIMENSION(klon,klev) :: exner,masse
119REAL, DIMENSION(klon,klev+1) :: masseb,q2old,q2neg
120
121!$OMP THREADPRIVATE(firstcall)       
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
140if (okiophys==1) then
141call iophys_ecrit('q2i',klev,'q2 debut my','m2/s2',q2(:,1:klev))
142call iophys_ecrit('kmi',klev,'Kz debut my','m/s2',km(:,1:klev))
143endif
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
204if (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))
218endif
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
340if (okiophys==1) then
341call iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev))
342call iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev))
343call iophys_ecrit('Km2',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))
344call iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev))
345endif
346#endif
347
348
349IF (iflag_pbl<20) then
350      ! For diagnostics only
351      RETURN
352
353ELSE
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
435if (okiophys==1) then
436call iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev))
437call iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev))
438call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
439call iophys_ecrit('q2neg',klev,'KTE non bornee','m2/s',q2neg(:,1:klev))
440call iophys_ecrit('alpha',klev,'alpha','',alpha(:,1:klev))
441call iophys_ecrit('sm',klev,'sm','',sm(:,1:klev))
442call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev))
443call iophys_ecrit('kmf',klev,'Kz final','m2/s',km(:,1:klev))
444call iophys_ecrit('knf',klev,'Kz final','m2/s',kn(:,1:klev))
445call iophys_ecrit('kqf',klev,'Kz final','m2/s',kq(:,1:klev))
446endif
447#endif
448
449ENDIF
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         stop'Cas nom prevu dans yamada4'
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
Note: See TracBrowser for help on using the repository browser.