source: LMDZ6/branches/Amaury_dev/libf/phylmd/yamada4.F90 @ 5218

Last change on this file since 5218 was 5117, checked in by abarral, 5 months ago

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