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

Last change on this file since 2828 was 2828, checked in by fhourdin, 7 years ago

Modifications de yamada4 issues de la these d'Etienne Vignon (sous

l'otpion new_yamada4=y) :
1/ schema numerique inspire de celui de MAR
2/ viscosite moleculaire ajoute pour les cas ou la diffusion

turbulente s'effondre.

3/ Desactivation de la modification Holtzlag Boville pour les

couches limites tres stables.

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