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

Last change on this file since 3774 was 3531, checked in by Laurent Fairhead, 5 years ago

Replaced STOP statements by a call to abort_physic in phylmd as per ticket #86
Still some work to be done in phylmd subdirectories

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