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

Last change on this file since 2598 was 2574, checked in by fhourdin, 8 years ago

Bug fixing

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