source: LMDZ5/trunk/libf/phylmd/yamada4.F90 @ 2722

Last change on this file since 2722 was 2721, checked in by fhourdin, 8 years ago

Modification de yamada4 pour corriger quelques erreurs et
autoriser de jouer avec les fonctions de stabilité.
Par defaut
new_yamada4=n donne les mêmes resultats numeriques qu'avant.
new_yamada4=y permet d'activer une nombre de Richardson critique
a partir du quel on ne fait plus décroitre les fonction de
stabilité Sm et Alpha de Mellor et Yamada.

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