source: LMDZ6/trunk/libf/phylmd/yamada4.F90 @ 3817

Last change on this file since 3817 was 3784, checked in by evignon, 4 years ago

correction v3783 pour convergence et compilation ancienne physique, Etienne aide par Ehouarn
  • 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
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.7 KB
Line 
1!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2
3SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
4    cd, tke, km, kn, kq, ustar, iflag_pbl, drgpro)
5
6  USE dimphy
7  USE ioipsl_getin_p_mod, ONLY : getin_p
8  USE phys_local_var_mod, only: tke_dissip
9 
10  IMPLICIT NONE
11  include "iniprint.h"
12  ! .......................................................................
13  ! ym#include "dimensions.h"
14  ! ym#include "dimphy.h"
15  ! ************************************************************************************************
16  !
17  ! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5
18  !
19  ! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model
20  !            Yamada T.
21  !            J. Atmos. Sci, 40, 91-106, 1983
22  !
23  !************************************************************************************************
24  ! Input :
25  !'======
26  ! ni: indice horizontal sur la grille de base, non restreinte
27  ! nsrf: type de surface
28  ! ngrid: nombre de mailles concern??es sur l'horizontal
29  ! dt : pas de temps
30  ! g  : g
31  ! rconst: constante de l'air sec
32  ! zlev : altitude a chaque niveau (interface inferieure de la couche
33  ! de meme indice)
34  ! zlay : altitude au centre de chaque couche
35  ! u,v : vitesse au centre de chaque couche
36  ! (en entree : la valeur au debut du pas de temps)
37  ! teta : temperature potentielle virtuelle au centre de chaque couche
38  ! (en entree : la valeur au debut du pas de temps)
39  ! cd : cdrag pour la quantit?? de mouvement
40  ! (en entree : la valeur au debut du pas de temps)
41  ! ustar: vitesse de friction calcul??e par une formule diagnostique
42  ! iflag_pbl: flag pour choisir des options du sch??ma de turbulence
43  !
44  !             iflag_pbl doit valoir entre 6 et 9
45  !             l=6, on prend  systematiquement une longueur d'equilibre
46  !             iflag_pbl=6 : MY 2.0
47  !             iflag_pbl=7 : MY 2.0.Fournier
48  !             iflag_pbl=8/9 : MY 2.5
49  !             iflag_pbl=8 with special obsolete treatments for convergence
50  !             with Cmpi5 NPv3.1 simulations
51  !             iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
52  !             iflag_pbl=12 = 11 with vertical diffusion off q2
53  !
54  !             2013/04/01 (FH hourdin@lmd.jussieu.fr)
55  !             Correction for very stable PBLs (iflag_pbl=10 and 11)
56  !             iflag_pbl=8 converges numerically with NPv3.1
57  !             iflag_pbl=11 -> the model starts with NP from start files created by ce0l
58  !                          -> the model can run with longer time-steps.
59  !             2016/11/30 (EV etienne.vignon@lmd.ipsl.fr)
60  !               On met tke (=q2/2) en entr??e plut??t que q2
61  !               On corrige l'update de la tke
62  !             2020/10/01 (EV)
63  !               On ajoute la dissipation de la TKE en diagnostique de sortie
64  !                 
65  ! Inpout/Output :
66  !==============
67  ! tke : tke au bas de chaque couche
68  ! (en entree : la valeur au debut du pas de temps)
69  ! (en sortie : la valeur a la fin du pas de temps)
70 
71  ! Outputs:
72  !==========
73  ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
74  ! couche)
75  ! (en sortie : la valeur a la fin du pas de temps)
76  ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
77  ! (en sortie : la valeur a la fin du pas de temps)
78  !
79  !.......................................................................
80
81  !=======================================================================
82  ! Declarations:
83  !=======================================================================
84
85
86  ! Inputs/Outputs
87  !----------------
88  REAL dt, g, rconst
89  REAL plev(klon, klev+1), temp(klon, klev)
90  REAL ustar(klon)
91  REAL kmin, qmin, pblhmin(klon), coriol(klon)
92  REAL zlev(klon, klev+1)
93  REAL zlay(klon, klev)
94  REAL u(klon, klev)
95  REAL v(klon, klev)
96  REAL teta(klon, klev)
97  REAL cd(klon)
98  REAL tke(klon, klev+1)
99  REAL unsdz(klon, klev)
100  REAL unsdzdec(klon, klev+1)
101  REAL kn(klon, klev+1)
102  REAL km(klon, klev+1)
103  INTEGER iflag_pbl, ngrid, nsrf
104  INTEGER ni(klon)
105
106!FC
107  REAL drgpro(klon,klev)
108  REAL winds(klon,klev)
109
110  ! Local
111  !-------
112
113  INCLUDE "clesphys.h"
114
115  REAL q2(klon, klev+1)
116  REAL kmpre(klon, klev+1), tmp2, qpre
117  REAL mpre(klon, klev+1)
118  REAL kq(klon, klev+1)
119  REAL ff(klon, klev+1), delta(klon, klev+1)
120  REAL aa(klon, klev+1), aa0, aa1
121  INTEGER nlay, nlev
122
123  LOGICAL,SAVE :: hboville=.TRUE.
124  REAL,SAVE :: viscom,viscoh
125  !$OMP THREADPRIVATE( hboville,viscom,viscoh)
126  INTEGER ig, jg, k
127  REAL ri, zrif, zalpha, zsm, zsn
128  REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
129  REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
130  REAL dtetadz(klon, klev+1)
131  REAL m2cstat, mcstat, kmcstat
132  REAL l(klon, klev+1)
133  REAL zz(klon, klev+1)
134  INTEGER iter
135  REAL dissip(klon,klev), tkeprov,tkeexp, shear(klon,klev), buoy(klon,klev)
136  REAL :: disseff
137
138  REAL,SAVE :: ric0,ric,rifc, b1, kap
139  !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
140  DATA b1, kap/16.6, 0.4/
141  REAL,SAVE :: seuilsm, seuilalpha
142  !$OMP THREADPRIVATE(seuilsm, seuilalpha)
143  REAL,SAVE :: lmixmin
144  !$OMP THREADPRIVATE(lmixmin)
145  LOGICAL, SAVE :: new_yamada4
146  INTEGER, SAVE :: yamada4_num
147  !$OMP THREADPRIVATE(new_yamada4,yamada4_num)
148  REAL, SAVE :: yun,ydeux
149  !$OMP THREADPRIVATE(yun,ydeux)
150
151  REAL frif, falpha, fsm
152  REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
153    lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
154  LOGICAL, SAVE :: firstcall = .TRUE.
155  !$OMP THREADPRIVATE(firstcall)
156
157  CHARACTER (len = 20) :: modname = 'yamada4'
158  CHARACTER (len = 80) :: abort_message
159
160
161
162  ! Fonctions utilis??es
163  !--------------------
164
165  frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
166  falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
167  fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
168 
169
170  IF (firstcall) THEN
171! Seuil dans le code de turbulence
172    new_yamada4=.false.
173    CALL getin_p('new_yamada4',new_yamada4)
174
175    IF (new_yamada4) THEN
176! Corrections et reglages issus du travail de these d'Etienne Vignon.
177       ric=0.143 ! qui donne des valeurs proches des seuils proposes
178                 ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
179                 !                    sm=1.1213 (1.12  dans Y83)
180       CALL getin_p('yamada4_ric',ric)
181       ric0=0.19489      ! ric=0.195 originalement, mais produisait sm<0
182       ric=min(ric,ric0) ! Au dela de ric0, sm devient n??gatif
183       rifc=frif(ric)
184       seuilsm=fsm(frif(ric))
185       seuilalpha=falpha(frif(ric))
186       yun=1.
187       ydeux=2.
188       hboville=.FALSE.
189       viscom=1.46E-5
190       viscoh=2.06E-5
191       !lmixmin=1.0E-3
192       lmixmin=0.
193       yamada4_num=5
194    ELSE
195       ric=0.195
196       rifc=0.191
197       seuilalpha=1.12
198       seuilsm=0.085
199       yun=2.
200       ydeux=1.
201       hboville=.TRUE.
202       viscom=0.
203       viscoh=0.
204       lmixmin=1.
205       yamada4_num=0
206    ENDIF
207
208    WRITE(lunout,*)'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
209    firstcall = .FALSE.
210    CALL getin_p('lmixmin',lmixmin)
211    CALL getin_p('yamada4_hboville',hboville)
212    CALL getin_p('yamada4_num',yamada4_num)
213  END IF
214
215
216
217!===============================================================================
218! Flags, tests et ??valuations de constantes
219!===============================================================================
220
221! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
222
223
224  IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN
225    abort_message='probleme de coherence dans appel a MY'
226    CALL abort_physic(modname,abort_message,1)
227  END IF
228
229
230  nlay = klev
231  nlev = klev + 1
232
233
234!========================================================================
235! Calcul des increments verticaux
236!=========================================================================
237
238 
239! Attention: zlev n'est pas declare a nlev
240  DO ig = 1, ngrid
241    zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
242  END DO
243
244
245  DO k = 1, nlay
246    DO ig = 1, ngrid
247      unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k))
248    END DO
249  END DO
250  DO ig = 1, ngrid
251    unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1))
252  END DO
253  DO k = 2, nlay
254    DO ig = 1, ngrid
255      unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1))
256    END DO
257  END DO
258  DO ig = 1, ngrid
259    unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
260  END DO
261
262!=========================================================================
263! Richardson number and stability functions
264!=========================================================================
265 
266! initialize arrays:
267
268  m2(1:ngrid, :) = 0.0
269  sm(1:ngrid, :) = 0.0
270  rif(1:ngrid, :) = 0.0
271
272!------------------------------------------------------------
273  DO k = 2, klev
274
275    DO ig = 1, ngrid
276      dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
277      m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
278        k-1))**2)/(dz(ig,k)*dz(ig,k))
279      dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
280      n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
281      ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
282      IF (ri<ric) THEN
283        rif(ig, k) = frif(ri)
284      ELSE
285        rif(ig, k) = rifc
286      END IF
287if (new_yamada4) then
288        alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha)
289        sm(ig, k) = max(fsm(rif(ig,k)),seuilsm)
290else
291      IF (rif(ig,k)<0.16) THEN
292        alpha(ig, k) = falpha(rif(ig,k))
293        sm(ig, k) = fsm(rif(ig,k))
294      ELSE
295        alpha(ig, k) = seuilalpha
296        sm(ig, k) = seuilsm
297      END IF
298
299end if
300      zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
301    END DO
302  END DO
303
304
305
306
307
308  !=======================================================================
309  !     DIFFERENT TYPES  DE SCHEMA de  YAMADA
310  !=======================================================================
311
312  ! On commence par calculer q2 a partir de la tke
313
314  IF (new_yamada4) THEN
315      DO k=1,klev+1
316         q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux
317      ENDDO
318  ELSE
319      DO k=1,klev+1
320         q2(1:ngrid,k)=tke(1:ngrid,k)
321      ENDDO
322  ENDIF
323
324! ====================================================================
325! Computing the mixing length
326! ====================================================================
327
328 
329  CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l)
330
331
332  !--------------
333  ! Yamada 2.0
334  !--------------
335  IF (iflag_pbl==6) THEN
336 
337    DO k = 2, klev
338      q2(1:ngrid, k) = l(1:ngrid, k)**2*zz(1:ngrid, k)
339    END DO
340
341
342  !------------------
343  ! Yamada 2.Fournier
344  !------------------
345
346  ELSE IF (iflag_pbl==7) THEN
347
348
349    ! Calcul de l,  km, au pas precedent
350    !....................................
351    DO k = 2, klev
352      DO ig = 1, ngrid
353        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
354        kmpre(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
355        mpre(ig, k) = sqrt(m2(ig,k))
356      END DO
357    END DO
358
359    DO k = 2, klev - 1
360      DO ig = 1, ngrid
361        m2cstat = max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1, 1.E-12)
362        mcstat = sqrt(m2cstat)
363
364     ! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3
365     !.........................................................................
366
367        IF (k==2) THEN
368          kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &
369            unsdz(ig,k-1)*cd(ig)*(sqrt(u(ig,3)**2+v(ig,3)**2)-mcstat/unsdzdec &
370            (ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k &
371            -1))
372        ELSE
373          kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &
374            unsdz(ig,k-1)*kmpre(ig,k-1)*mpre(ig,k-1))/ &
375            (unsdz(ig,k)+unsdz(ig,k-1))
376        END IF
377
378        tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
379        q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
380
381      END DO
382    END DO
383
384
385    ! ------------------------
386    ! Yamada 2.5 a la Didi
387    !-------------------------
388
389  ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
390
391    ! Calcul de l, km, au pas precedent
392    !....................................
393    DO k = 2, klev
394      DO ig = 1, ngrid
395        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
396        IF (delta(ig,k)<1.E-20) THEN
397          delta(ig, k) = 1.E-20
398        END IF
399        km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
400        aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
401        aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
402        aa(ig, k) = aa1*dt/(delta(ig,k)*l(ig,k))
403        qpre = sqrt(q2(ig,k))
404        IF (aa(ig,k)>0.) THEN
405          q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2
406        ELSE
407          q2(ig, k) = (qpre/(1.-aa(ig,k)*qpre))**2
408        END IF
409        ! else ! iflag_pbl=9
410        ! if (aa(ig,k)*qpre.gt.0.9) then
411        ! q2(ig,k)=(qpre*10.)**2
412        ! else
413        ! q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
414        ! endif
415        ! endif
416        q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
417      END DO
418    END DO
419
420  ELSE IF (iflag_pbl>=10) THEN
421
422    shear(:,:)=0.
423    buoy(:,:)=0.
424    dissip(:,:)=0.
425    km(:,:)=0.
426       
427    IF (yamada4_num>=1) THEN
428 
429    DO k = 2, klev - 1
430      DO ig=1,ngrid
431      q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
432      km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
433      shear(ig,k)=km(ig, k)*m2(ig, k)
434      buoy(ig,k)=km(ig, k)*m2(ig, k)*(-1.*rif(ig,k))
435!      dissip(ig,k)=min(max(((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)),1.E-12),1.E4)
436      dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k))
437     ENDDO
438    ENDDO
439   
440    IF (yamada4_num==1) THEN ! Schema du MAR tel quel
441       DO k = 2, klev - 1
442         DO ig=1,ngrid
443         tkeprov=q2(ig,k)/ydeux
444         tkeprov= tkeprov*                           &
445           &  (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ &
446           &  (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
447         q2(ig,k)=tkeprov*ydeux
448        ENDDO
449       ENDDO
450    ELSE IF (yamada4_num==2) THEN ! version modifiee avec integration exacte pour la dissipation
451       DO k = 2, klev - 1
452         DO ig=1,ngrid
453         tkeprov=q2(ig,k)/ydeux
454         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
455         tkeprov = tkeprov/(1.+dt*disseff/(2.*tkeprov))**2
456         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
457         q2(ig,k)=tkeprov*ydeux
458         ! En cas stable, on traite la flotabilite comme la
459         ! dissipation, en supposant que buoy/q2^3 est constant.
460         ! Puis on prend la solution exacte
461        ENDDO
462       ENDDO
463    ELSE IF (yamada4_num==3) THEN ! version modifiee avec integration exacte pour la dissipation
464       DO k = 2, klev - 1
465         DO ig=1,ngrid
466         tkeprov=q2(ig,k)/ydeux
467         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
468         tkeprov=tkeprov*exp(-dt*disseff/tkeprov)
469         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
470         q2(ig,k)=tkeprov*ydeux
471         ! En cas stable, on traite la flotabilite comme la
472         ! dissipation, en supposant que buoy/q2^3 est constant.
473         ! Puis on prend la solution exacte
474        ENDDO
475       ENDDO
476    ELSE IF (yamada4_num==4) THEN ! version modifiee avec integration exacte pour la dissipation
477       DO k = 2, klev - 1
478         DO ig=1,ngrid
479         tkeprov=q2(ig,k)/ydeux
480         tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
481         tkeprov= tkeprov*                           &
482           &  tkeprov/ &
483           &  (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
484         q2(ig,k)=tkeprov*ydeux
485         ! En cas stable, on traite la flotabilite comme la
486         ! dissipation, en supposant que buoy/q2^3 est constant.
487         ! Puis on prend la solution exacte
488        ENDDO
489       ENDDO
490       
491      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
492      !! Attention, yamada4_num=5 est inexacte car néglige les termes de flottabilité
493      !! en conditions instables
494      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
495    ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation
496       DO k = 2, klev - 1
497         DO ig=1,ngrid
498         tkeprov=q2(ig,k)/ydeux
499
500!             if(ifl_pbltree .eq. 0) then
501!         disseff=dissip(ig,k)-min(0.,buoy(ig,k))
502!         tkeexp=exp(-dt*disseff/tkeprov)
503!         tkeprov= shear(ig,k)*tkeprov/disseff*(1.-tkeexp)+tkeprov*tkeexp
504!           else
505!FC on ajoute la dissipation due aux arbres
506         disseff=dissip(ig,k)-min(0.,buoy(ig,k)) + drgpro(ig,k)*tkeprov
507         tkeexp=exp(-dt*disseff/tkeprov)
508! on prend en compte la tke cree par les arbres
509         winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
510         tkeprov= (shear(ig,k)+ &
511          & drgpro(ig,k)*(winds(ig,k))**3)*tkeprov/disseff*(1.-tkeexp)+tkeprov*tkeexp
512!               endif
513
514         q2(ig,k)=tkeprov*ydeux
515
516         ! En cas stable, on traite la flotabilite comme la
517         ! dissipation, en supposant que buoy/q2^3 est constant.
518         ! Puis on prend la solution exacte
519        ENDDO
520       ENDDO
521    ELSE IF (yamada4_num==6) THEN ! version modifiee avec integration exacte pour la dissipation
522       DO k = 2, klev - 1
523         DO ig=1,ngrid
524         !tkeprov=q2(ig,k)/ydeux
525         !tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt
526         !disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k))
527         !tkeexp=exp(-dt*disseff/tkeprov)
528         !tkeprov= tkeprov*tkeexp
529         !q2(ig,k)=tkeprov*ydeux
530         ! En cas stable, on traite la flotabilite comme la
531         ! dissipation, en supposant que dissipeff/TKE est constant.
532         ! Puis on prend la solution exacte
533         !
534         ! With drag and dissipation from high vegetation (EV & FC, 05/10/2020)
535         winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
536         tkeprov=q2(ig,k)/ydeux
537         tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3,0.)*dt
538         disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3) + drgpro(ig,k)*tkeprov
539         tkeexp=exp(-dt*disseff/tkeprov)
540         tkeprov= tkeprov*tkeexp
541         q2(ig,k)=tkeprov*ydeux
542         
543        ENDDO
544       ENDDO
545    ENDIF
546
547    DO k = 2, klev - 1
548      DO ig=1,ngrid
549      q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
550      ENDDO
551    ENDDO
552
553   ELSE
554
555    DO k = 2, klev - 1
556      km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k)
557      q2(1:ngrid, k) = q2(1:ngrid, k) + ydeux*dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
558!     q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
559      q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
560       q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1))
561!     q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
562      q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
563    END DO
564
565  ENDIF
566
567  ELSE
568     abort_message='Cas nom prevu dans yamada4'
569     CALL abort_physic(modname,abort_message,1)
570
571  END IF ! Fin du cas 8
572
573
574  ! ====================================================================
575  ! Calcul des coefficients de melange
576  ! ====================================================================
577
578  DO k = 2, klev
579    DO ig = 1, ngrid
580      zq = sqrt(q2(ig,k))
581      km(ig, k) = l(ig, k)*zq*sm(ig, k)     ! For momentum
582      kn(ig, k) = km(ig, k)*alpha(ig, k)    ! For scalars
583      kq(ig, k) = l(ig, k)*zq*0.2           ! For TKE
584    END DO
585  END DO
586
587
588  !====================================================================
589  ! Transport diffusif vertical de la TKE par la TKE
590  !====================================================================
591
592
593    ! initialize near-surface and top-layer mixing coefficients
594    !...........................................................
595
596  kq(1:ngrid, 1) = kq(1:ngrid, 2)    ! constant (ie no gradient) near the surface
597  kq(1:ngrid, klev+1) = 0            ! zero at the top
598
599    ! Transport diffusif vertical de la TKE.
600    !.......................................
601
602  IF (iflag_pbl>=12) THEN
603    q2(1:ngrid, 1) = q2(1:ngrid, 2)
604    CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
605  END IF
606
607
608  !====================================================================
609  ! Traitement particulier pour les cas tres stables, introduction d'une
610  ! longueur de m??lange minimale
611  !====================================================================
612  !
613  ! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model
614  !            Holtslag A.A.M. and Boville B.A.
615  !            J. Clim., 6, 1825-1842, 1993
616
617
618 IF (hboville) THEN
619
620
621  IF (prt_level>1) THEN
622    WRITE(lunout,*) 'YAMADA4 0'
623  END IF
624
625  DO ig = 1, ngrid
626    coriol(ig) = 1.E-4
627    pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
628  END DO
629
630  IF (1==1) THEN
631    IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
632
633      DO k = 2, klev
634        DO ig = 1, ngrid
635          IF (teta(ig,2)>teta(ig,1)) THEN
636            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
637            kmin = kap*zlev(ig, k)*qmin
638          ELSE
639            kmin = -1. ! kmin n'est utilise que pour les SL stables.
640          END IF
641          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
642
643            kn(ig, k) = kmin
644            km(ig, k) = kmin
645            kq(ig, k) = kmin
646
647 ! la longueur de melange est suposee etre l= kap z
648 ! K=l q Sm d'ou q2=(K/l Sm)**2
649
650            q2(ig, k) = (qmin/sm(ig,k))**2
651          END IF
652        END DO
653      END DO
654
655    ELSE
656      DO k = 2, klev
657        DO ig = 1, ngrid
658          IF (teta(ig,2)>teta(ig,1)) THEN
659            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
660            kmin = kap*zlev(ig, k)*qmin
661          ELSE
662            kmin = -1. ! kmin n'est utilise que pour les SL stables.
663          END IF
664          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
665            kn(ig, k) = kmin
666            km(ig, k) = kmin
667            kq(ig, k) = kmin
668 ! la longueur de melange est suposee etre l= kap z
669 ! K=l q Sm d'ou q2=(K/l Sm)**2
670            sm(ig, k) = 1.
671            alpha(ig, k) = 1.
672            q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
673            zq = sqrt(q2(ig,k))
674            km(ig, k) = l(ig, k)*zq*sm(ig, k)
675            kn(ig, k) = km(ig, k)*alpha(ig, k)
676            kq(ig, k) = l(ig, k)*zq*0.2
677          END IF
678        END DO
679      END DO
680    END IF
681
682  END IF
683
684 END IF ! hboville
685
686! Ajout d'une viscosite moleculaire
687   km(1:ngrid,2:klev)=km(1:ngrid,2:klev)+viscom
688   kn(1:ngrid,2:klev)=kn(1:ngrid,2:klev)+viscoh
689   kq(1:ngrid,2:klev)=kq(1:ngrid,2:klev)+viscoh
690
691  IF (prt_level>1) THEN
692    WRITE(lunout,*)'YAMADA4 1'
693  END IF !(prt_level>1) THEN
694
695
696 !======================================================
697 ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
698 !======================================================
699 !
700 ! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer
701 !            Abdella K and McFarlane N
702 !            J. Atmos. Sci., 54, 1850-1867, 1997
703
704  ! Diagnostique pour stokage
705  !..........................
706
707  IF (1==0) THEN
708    rino = rif
709    smyam(1:ngrid, 1) = 0.
710    styam(1:ngrid, 1) = 0.
711    lyam(1:ngrid, 1) = 0.
712    knyam(1:ngrid, 1) = 0.
713    w2yam(1:ngrid, 1) = 0.
714    t2yam(1:ngrid, 1) = 0.
715
716    smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
717    styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
718    lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
719    knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
720
721
722  ! Calcul de w'2 et T'2
723  !.......................
724
725    w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
726      lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
727      sqrt(q2(1:ngrid,2:klev))
728
729    t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* &
730      dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* &
731      lyam(1:ngrid, 2:klev)
732  END IF
733
734
735
736!============================================================================
737! Mise a jour de la tke
738!============================================================================
739
740  IF (new_yamada4) THEN
741     DO k=1,klev+1
742        tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux
743     ENDDO
744  ELSE
745     DO k=1,klev+1
746        tke(1:ngrid,k)=q2(1:ngrid,k)
747     ENDDO
748  ENDIF
749
750
751!============================================================================
752! Diagnostique de la dissipation
753!============================================================================
754
755! Diagnostics
756 tke_dissip(1:ngrid,:,nsrf)=0.
757! DO k=2,klev
758!    DO ig=1,ngrid
759!       jg=ni(ig)
760!       tke_dissip(jg,k,nsrf)=dissip(ig,k)
761!    ENDDO
762! ENDDO
763 
764!=============================================================================
765
766  RETURN
767
768
769END SUBROUTINE yamada4
770
771!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
772
773!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
774SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
775
776  USE dimphy
777  IMPLICIT NONE
778 
779!    vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE
780!             avec un schema implicite en temps avec
781!             inversion d'un syst??me tridiagonal
782!
783!     Reference: Description of the interface with the surface and
784!                the computation of the turbulet diffusion in LMDZ
785!                Technical note on LMDZ
786!                Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y
787!
788!============================================================================
789! Declarations
790!============================================================================
791
792  REAL plev(klon, klev+1)
793  REAL temp(klon, klev)
794  REAL timestep
795  REAL gravity, rconst
796  REAL kstar(klon, klev+1), zz
797  REAL kmy(klon, klev+1)
798  REAL q2(klon, klev+1)
799  REAL deltap(klon, klev+1)
800  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
801  INTEGER ngrid
802
803  INTEGER i, k
804
805
806!=========================================================================
807! Calcul
808!=========================================================================
809
810  DO k = 1, klev
811    DO i = 1, ngrid
812      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
813      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
814        (plev(i,k)-plev(i,k+1))*timestep
815    END DO
816  END DO
817
818  DO k = 2, klev
819    DO i = 1, ngrid
820      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
821    END DO
822  END DO
823  DO i = 1, ngrid
824    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
825    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
826    denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev)
827    alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1)
828    beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1)
829  END DO
830
831  DO k = klev, 2, -1
832    DO i = 1, ngrid
833      denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
834        kstar(i, k-1)
835      alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k)
836      beta(i, k) = kstar(i, k-1)/denom(i, k)
837    END DO
838  END DO
839
840  ! Si on recalcule q2(1)
841  !.......................
842  IF (1==0) THEN
843    DO i = 1, ngrid
844      denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)
845      q2(i, 1) = (q2(i,1)*deltap(i,1)+kstar(i,1)*alpha(i,2))/denom(i, 1)
846    END DO
847  END IF
848
849
850  DO k = 2, klev + 1
851    DO i = 1, ngrid
852      q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1)
853    END DO
854  END DO
855
856  RETURN
857END SUBROUTINE vdif_q2
858!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
859
860
861
862!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
863 SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
864 
865   USE dimphy
866  IMPLICIT NONE
867
868! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE
869!           avec un schema explicite en temps
870
871
872!====================================================
873! Declarations
874!====================================================
875
876  REAL plev(klon, klev+1)
877  REAL temp(klon, klev)
878  REAL timestep
879  REAL gravity, rconst
880  REAL kstar(klon, klev+1), zz
881  REAL kmy(klon, klev+1)
882  REAL q2(klon, klev+1)
883  REAL deltap(klon, klev+1)
884  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
885  INTEGER ngrid
886  INTEGER i, k
887
888
889!==================================================
890! Calcul
891!==================================================
892
893  DO k = 1, klev
894    DO i = 1, ngrid
895      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
896      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
897        (plev(i,k)-plev(i,k+1))*timestep
898    END DO
899  END DO
900
901  DO k = 2, klev
902    DO i = 1, ngrid
903      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
904    END DO
905  END DO
906  DO i = 1, ngrid
907    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
908    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
909  END DO
910
911  DO k = klev, 2, -1
912    DO i = 1, ngrid
913      q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
914        k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
915    END DO
916  END DO
917
918  DO i = 1, ngrid
919    q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
920    q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
921      klev)))/deltap(i, klev+1)
922  END DO
923
924  RETURN
925END SUBROUTINE vdif_q2e
926
927!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
928
929
930!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
931
932SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix)
933
934
935
936  USE dimphy
937  USE phys_state_var_mod, only: zstd, zsig, zmea
938  USE phys_local_var_mod, only: l_mixmin, l_mix
939
940 ! zstd: ecart type de la'altitud e sous-maille
941 ! zmea: altitude moyenne sous maille
942 ! zsig: pente moyenne de le maille
943
944  USE geometry_mod, only: cell_area
945  ! aire_cell: aire de la maille
946
947  IMPLICIT NONE
948!*************************************************************************
949! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence
950! avec la formule de Blackadar
951! Calcul d'un  minimum en fonction de l'orographie sous-maille:
952! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange
953! induit par les circulations meso et submeso ??chelles.
954!
955! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere
956!               Blackadar A.K.
957!               J. Geophys. Res., 64, No 8, 1962
958!
959!             * An evaluation of neutral and convective planetary boundary-layer parametrisations relative
960!               to large eddy simulations
961!               Ayotte K et al
962!               Boundary Layer Meteorology, 79, 131-175, 1996
963!
964!
965!             * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts
966!               Van de Wiel B.J.H et al
967!               Boundary-Lay Meteorol, 128, 103-166, 2008
968!
969!
970! Histoire:
971!----------
972! * premi??re r??daction, Etienne et Frederic, 09/06/2016
973!
974! ***********************************************************************
975
976!==================================================================
977! Declarations
978!==================================================================
979
980! Inputs
981!-------
982 INTEGER            ni(klon)           ! indice sur la grille original (non restreinte)
983 INTEGER            nsrf               ! Type de surface
984 INTEGER            ngrid              ! Nombre de points concern??s sur l'horizontal
985 INTEGER            iflag_pbl          ! Choix du sch??ma de turbulence
986 REAL            pbl_lmixmin_alpha  ! on active ou non le calcul de la longueur de melange minimum
987 REAL               lmixmin            ! Minimum absolu de la longueur de m??lange
988 REAL               zlay(klon, klev)   ! altitude du centre de la couche
989 REAL               zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche
990 REAL               u(klon, klev)      ! vitesse du vent zonal
991 REAL               v(klon, klev)      ! vitesse du vent meridional
992 REAL               q2(klon, klev+1)   ! energie cin??tique turbulente
993 REAL               n2(klon, klev+1)   ! frequence de Brunt-Vaisala
994
995!In/out
996!-------
997
998  LOGICAL, SAVE :: firstcall = .TRUE.
999  !$OMP THREADPRIVATE(firstcall)
1000
1001! Outputs
1002!---------
1003
1004 REAL               lmix(klon, klev+1)    ! Longueur de melange 
1005
1006
1007! Local
1008!-------
1009 
1010 INTEGER  ig,jg, k
1011 REAL     h_oro(klon)
1012 REAL     hlim(klon)
1013 REAL, SAVE :: kap=0.4,kapb=0.4
1014  !$OMP THREADPRIVATE(kap,kapb)
1015 REAL zq
1016 REAL sq(klon), sqz(klon)
1017 REAL, ALLOCATABLE, SAVE :: l0(:)
1018  !$OMP THREADPRIVATE(l0)
1019 REAL fl, zzz, zl0, zq2, zn2
1020 REAL famorti, zzzz, zh_oro, zhlim
1021 REAL l1(klon, klev+1), l2(klon,klev+1)
1022 REAL winds(klon, klev)
1023 REAL xcell
1024 REAL zstdslope(klon) 
1025 REAL lmax
1026 REAL l2strat, l2neutre, extent 
1027 REAL l2limit(klon)
1028!===============================================================
1029! Fonctions utiles
1030!===============================================================
1031
1032! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996
1033!..........................................................................
1034
1035 fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &
1036    k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &
1037    max(n2(ig,k),1.E-10))), 1.E-5)
1038 
1039! Fonction d'amortissement de la turbulence au dessus de la montagne
1040! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code
1041!.....................................................................
1042
1043 famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1.   
1044
1045  IF (ngrid==0) RETURN
1046
1047  IF (firstcall) THEN
1048    ALLOCATE (l0(klon))
1049    firstcall = .FALSE.
1050  END IF
1051
1052
1053!=====================================================================
1054!         CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1
1055!=====================================================================
1056
1057  l1(1:ngrid,:)=0.
1058  IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
1059
1060   
1061    ! Iterative computation of l0
1062    ! This version is kept for iflag_pbl only for convergence
1063    ! with NPv3.1 Cmip5 simulations
1064    !...................................................................
1065
1066    DO ig = 1, ngrid
1067      sq(ig) = 1.E-10
1068      sqz(ig) = 1.E-10
1069    END DO
1070    DO k = 2, klev - 1
1071      DO ig = 1, ngrid
1072        zq = sqrt(q2(ig,k))
1073        sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
1074        sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
1075      END DO
1076    END DO
1077    DO ig = 1, ngrid
1078      l0(ig) = 0.2*sqz(ig)/sq(ig)
1079    END DO
1080    DO k = 2, klev
1081      DO ig = 1, ngrid
1082        l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
1083      END DO
1084    END DO
1085
1086  ELSE
1087
1088   
1089    ! In all other case, the assymptotic mixing length l0 is imposed (150m)
1090    !......................................................................
1091
1092    l0(1:ngrid) = 150.
1093    DO k = 2, klev
1094      DO ig = 1, ngrid
1095        l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
1096      END DO
1097    END DO
1098
1099  END IF
1100
1101!=================================================================================
1102!  CALCUL d'une longueur de melange en fonctions de la topographie sous maille: l2
1103! si plb_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les
1104! glacier, la glace de mer et les oc??ans)
1105!=================================================================================
1106
1107   l2(1:ngrid,:)=0.0
1108   l_mixmin(1:ngrid,:,nsrf)=0.
1109   l_mix(1:ngrid,:,nsrf)=0.
1110
1111   IF (nsrf .EQ. 1) THEN
1112
1113! coefficients
1114!--------------
1115
1116     extent=2.                                                         ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. 
1117     lmax=150.                                                         ! Longueur de m??lange max dans l'absolu
1118
1119! calculs
1120!---------
1121
1122     DO ig=1,ngrid
1123
1124      ! On calcule la hauteur du relief
1125      !.................................
1126      ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille
1127      ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon)
1128      ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief
1129      ! (en gros, une maille de taille xcell avec une pente constante zstdslope)
1130      jg=ni(ig)
1131!     IF (zsig(jg) .EQ. 0.) THEN
1132!          zstdslope(ig)=0.         
1133!     ELSE
1134!     xcell=sqrt(cell_area(jg))
1135!     zstdslope(ig)=max((xcell*zsig(jg)-zmea(jg))**3 /(3.*zsig(jg)),0.)
1136!     zstdslope(ig)=sqrt(zstdslope(ig))
1137!     END IF
1138     
1139!     h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.)   ! Hauteur du relief
1140      h_oro(ig)=zstd(jg)
1141      hlim(ig)=extent*h_oro(ig)     
1142     ENDDO
1143
1144     l2limit(1:ngrid)=0.
1145
1146     DO k=2,klev
1147        DO ig=1,ngrid
1148           winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
1149           IF (zlev(ig,k) .LE. h_oro(ig)) THEN  ! sous l'orographie
1150              l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10))  ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008)
1151              l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig))            ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h
1152              l2neutre=MIN(l2neutre,lmax)                                             ! On majore par lmax
1153              l2limit(ig)=MIN(l2neutre,l2strat)                                       ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e)
1154              l2(ig,k)=l2limit(ig)
1155                                     
1156           ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles
1157
1158      ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes
1159      ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z)
1160      ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim
1161              l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig))
1162           ELSE                                                                    ! Au dessus de extent*h, on prend l2=l0
1163              l2(ig,k)=0.
1164           END IF
1165        ENDDO
1166     ENDDO
1167   ENDIF                                                                        ! pbl_lmixmin_alpha
1168
1169!==================================================================================
1170! On prend le max entre la longueur de melange de blackadar et celle calcul??e
1171! en fonction de la topographie
1172!===================================================================================
1173
1174
1175 DO k=1,klev+1
1176    DO ig=1,ngrid
1177       lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
1178   ENDDO
1179 ENDDO
1180
1181! Diagnostics
1182
1183 DO k=2,klev
1184    DO ig=1,ngrid
1185       jg=ni(ig)
1186       l_mix(jg,k,nsrf)=lmix(ig,k)
1187       l_mixmin(jg,k,nsrf)=l2(ig,k)
1188    ENDDO
1189 ENDDO
1190 DO ig=1,ngrid
1191    jg=ni(ig)
1192    l_mix(jg,1,nsrf)=hlim(ig)
1193 ENDDO
1194
1195
1196
1197END SUBROUTINE mixinglength
Note: See TracBrowser for help on using the repository browser.