Sources/phylmd/yamada4.f

module yamada4_m

  IMPLICIT NONE

  real, parameter:: kap = 0.4
  private
  public yamada4

contains

  SUBROUTINE yamada4(ngrid, dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, &
       kq, ustar, iflag_pbl)

    ! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36

    USE dimphy, ONLY : klev, klon

    integer ngrid
    REAL, intent(in):: dt ! pas de temps
    real, intent(in):: g

    REAL zlev(klon, klev+1)
    ! altitude à chaque niveau (interface inférieure de la couche de
    ! même indice)

    REAL zlay(klon, klev) ! altitude au centre de chaque couche

    REAL u(klon, klev), v(klon, klev)
    ! vitesse au centre de chaque couche (en entrée : la valeur au
    ! début du pas de temps)

    REAL teta(klon, klev)
    ! température potentielle au centre de chaque couche (en entrée :
    ! la valeur au début du pas de temps)

    REAL cd(klon) ! cdrag (en entrée : la valeur au début du pas de temps)

    REAL, intent(inout):: q2(klon, klev+1)
    ! $q^2$ au bas de chaque couche 
    ! En entrée : la valeur au début du pas de temps ; en sortie : la
    ! valeur à la fin du pas de temps.

    REAL km(klon, klev+1)
    ! diffusivité turbulente de quantité de mouvement (au bas de
    ! chaque couche) (en sortie : la valeur à la fin du pas de temps)

    REAL kn(klon, klev+1)
    ! diffusivité turbulente des scalaires (au bas de chaque couche)
    ! (en sortie : la valeur à la fin du pas de temps)

    REAL kq(klon, klev+1)
    real ustar(klon)

    integer iflag_pbl
    ! iflag_pbl doit valoir entre 6 et 9
    ! l = 6, on prend systématiquement une longueur d'équilibre
    ! iflag_pbl = 6 : MY 2.0
    ! iflag_pbl = 7 : MY 2.0.Fournier
    ! iflag_pbl = 8 : MY 2.5
    ! iflag_pbl = 9 : un test ?

    ! Local:

    real kmin, qmin, pblhmin(klon), coriol(klon)
    real qpre
    REAL unsdz(klon, klev)
    REAL unsdzdec(klon, klev+1)
    REAL kmpre(klon, klev+1), tmp2
    REAL mpre(klon, klev+1)
    real delta(klon, klev+1)
    real aa(klon, klev+1), aa0, aa1
    integer, PARAMETER:: nlay = klev
    integer, PARAMETER:: nlev = klev+1
    logical:: first = .true.
    integer:: ipas = 0
    integer ig, k
    real ri
    real rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
    real m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
    real dtetadz(klon, klev+1)
    real m2cstat, mcstat, kmcstat
    real l(klon, klev+1)
    real, save:: l0(klon)
    real sq(klon), sqz(klon), zz(klon, klev+1)
    integer iter
    real:: ric = 0.195, rifc = 0.191, b1 = 16.6, kap = 0.4
    real rino(klon, klev+1), smyam(klon, klev), styam(klon, klev)
    real lyam(klon, klev), knyam(klon, klev)

    !-----------------------------------------------------------------------

    if (.not. (iflag_pbl >= 6 .and. iflag_pbl <= 9)) then
       print *, 'probleme de coherence dans appel a MY'
       stop 1
    endif

    ipas = ipas+1

    ! les increments verticaux
    DO ig = 1, ngrid
       ! alerte: zlev n'est pas declare a nlev
       zlev(ig, nlev) = zlay(ig, nlay) +(zlay(ig, nlay) - zlev(ig, nlev-1))
    ENDDO

    DO k = 1, nlay
       DO ig = 1, ngrid
          unsdz(ig, k) = 1.E+0/(zlev(ig, k+1)-zlev(ig, k))
       ENDDO
    ENDDO
    DO ig = 1, ngrid
       unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))
    ENDDO
    DO k = 2, nlay
       DO ig = 1, ngrid
          unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))
       ENDDO
    ENDDO
    DO ig = 1, ngrid
       unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig, nlay+1)-zlay(ig, nlay))
    ENDDO

    do k = 2, klev
       do ig = 1, ngrid
          dz(ig, k) = zlay(ig, k)-zlay(ig, k-1)
          m2(ig, k) = ((u(ig, k)-u(ig, k-1))**2+(v(ig, k)-v(ig, k-1))**2) &
               /(dz(ig, k)*dz(ig, k))
          dtetadz(ig, k) = (teta(ig, k)-teta(ig, k-1))/dz(ig, k)
          n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig, k-1)+teta(ig, k))
          ri = n2(ig, k)/max(m2(ig, k), 1.e-10)
          if (ri.lt.ric) then
             rif(ig, k) = frif(ri)
          else
             rif(ig, k) = rifc
          endif
          if(rif(ig, k).lt.0.16) then
             alpha(ig, k) = falpha(rif(ig, k))
             sm(ig, k) = fsm(rif(ig, k))
          else
             alpha(ig, k) = 1.12
             sm(ig, k) = 0.085
          endif
          zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig, k))*sm(ig, k)
       enddo
    enddo

    ! Au premier appel, on détermine l et q2 de façon itérative.
    ! Itération pour déterminer la longueur de mélange

    if (first .or. iflag_pbl == 6) then
       do ig = 1, ngrid
          l0(ig) = 10.
       enddo
       do k = 2, klev-1
          do ig = 1, ngrid
             l(ig, k) = l0(ig) * kap * zlev(ig, k) &
                  / (kap * zlev(ig, k) + l0(ig))
          enddo
       enddo

       do iter = 1, 10
          do ig = 1, ngrid
             sq(ig) = 1e-10
             sqz(ig) = 1e-10
          enddo
          do k = 2, klev-1
             do ig = 1, ngrid
                q2(ig, k) = l(ig, k)**2 * zz(ig, k)
                l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
                zq = sqrt(q2(ig, k))
                sqz(ig) = sqz(ig) + zq * zlev(ig, k) &
                     * (zlay(ig, k) - zlay(ig, k-1))
                sq(ig) = sq(ig) + zq * (zlay(ig, k) - zlay(ig, k-1))
             enddo
          enddo
          do ig = 1, ngrid
             l0(ig) = 0.2 * sqz(ig) / sq(ig)
          enddo
       enddo
    endif

    ! Calcul de la longueur de melange.

    ! Mise a jour de l0
    do ig = 1, ngrid
       sq(ig) = 1.e-10
       sqz(ig) = 1.e-10
    enddo
    do k = 2, klev-1
       do ig = 1, ngrid
          zq = sqrt(q2(ig, k))
          sqz(ig) = sqz(ig)+zq*zlev(ig, k)*(zlay(ig, k)-zlay(ig, k-1))
          sq(ig) = sq(ig)+zq*(zlay(ig, k)-zlay(ig, k-1))
       enddo
    enddo
    do ig = 1, ngrid
       l0(ig) = 0.2*sqz(ig)/sq(ig)
    enddo
    ! calcul de l(z)
    do k = 2, klev
       do ig = 1, ngrid
          l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
          if(first) then
             q2(ig, k) = l(ig, k)**2 * zz(ig, k)
          endif
       enddo
    enddo

    ! Yamada 2.0
    if (iflag_pbl == 6) then
       do k = 2, klev
          do ig = 1, ngrid
             q2(ig, k) = l(ig, k)**2 * zz(ig, k)
          enddo
       enddo
    else if (iflag_pbl == 7) then
       ! Yamada 2.Fournier

       ! Calcul de l, km, au pas precedent
       do k = 2, klev
          do ig = 1, ngrid
             delta(ig, k) = q2(ig, k) / (l(ig, k)**2 * sm(ig, k))
             kmpre(ig, k) = l(ig, k) * sqrt(q2(ig, k)) * sm(ig, k)
             mpre(ig, k) = sqrt(m2(ig, k))
          enddo
       enddo

       do k = 2, klev-1
          do ig = 1, ngrid
             m2cstat = max(alpha(ig, k)*n2(ig, k)+delta(ig, k)/b1, 1.e-12)
             mcstat = sqrt(m2cstat)

             ! puis on ecrit la valeur de q qui annule l'equation de m
             ! supposee en q3

             IF (k == 2) THEN
                kmcstat = 1.E+0 / mcstat &
                     *(unsdz(ig, k)*kmpre(ig, k+1) &
                     *mpre(ig, k+1) &
                     +unsdz(ig, k-1) &
                     *cd(ig) &
                     *(sqrt(u(ig, 3)**2+v(ig, 3)**2) &
                     -mcstat/unsdzdec(ig, k) &
                     -mpre(ig, k+1)/unsdzdec(ig, k+1))**2) &
                     /(unsdz(ig, k)+unsdz(ig, k-1))
             ELSE
                kmcstat = 1.E+0 / mcstat &
                     *(unsdz(ig, k)*kmpre(ig, k+1) &
                     *mpre(ig, k+1) &
                     +unsdz(ig, k-1)*kmpre(ig, k-1) &
                     *mpre(ig, k-1)) &
                     /(unsdz(ig, k)+unsdz(ig, k-1))
             ENDIF
             tmp2 = kmcstat / (sm(ig, k) / q2(ig, k)) /l(ig, k)
             q2(ig, k) = max(tmp2, 1.e-12)**(2./3.)
          enddo
       enddo
    else if (iflag_pbl >= 8) then
       ! Yamada 2.5 a la Didi

       ! Calcul de l, km, au pas precedent
       do k = 2, klev
          do ig = 1, ngrid
             delta(ig, k) = q2(ig, k)/(l(ig, k)**2*sm(ig, k))
             if (delta(ig, k).lt.1.e-20) then
                delta(ig, k) = 1.e-20
             endif
             km(ig, k) = l(ig, k)*sqrt(q2(ig, k))*sm(ig, k)
             aa0 = (m2(ig, k)-alpha(ig, k)*n2(ig, k)-delta(ig, k)/b1)
             aa1 = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1)
             aa(ig, k) = aa1*dt/(delta(ig, k)*l(ig, k))
             qpre = sqrt(q2(ig, k))
             if (iflag_pbl == 8) then
                if (aa(ig, k).gt.0.) then
                   q2(ig, k) = (qpre+aa(ig, k)*qpre*qpre)**2
                else
                   q2(ig, k) = (qpre/(1.-aa(ig, k)*qpre))**2
                endif
             else
                ! iflag_pbl = 9
                if (aa(ig, k)*qpre.gt.0.9) then
                   q2(ig, k) = (qpre*10.)**2
                else
                   q2(ig, k) = (qpre/(1.-aa(ig, k)*qpre))**2
                endif
             endif
             q2(ig, k) = min(max(q2(ig, k), 1.e-10), 1.e4)
          enddo
       enddo
    endif

    ! Calcul des coefficients de mélange
    do k = 2, klev
       do ig = 1, ngrid
          zq = sqrt(q2(ig, k))
          km(ig, k) = l(ig, k)*zq*sm(ig, k)
          kn(ig, k) = km(ig, k)*alpha(ig, k)
          kq(ig, k) = l(ig, k)*zq*0.2
       enddo
    enddo

    ! Traitement des cas noctrunes avec l'introduction d'une longueur
    ! minilale.

    ! Traitement particulier pour les cas tres stables.
    ! D'apres Holtslag Boville.

    do ig = 1, ngrid
       coriol(ig) = 1.e-4
       pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546e-5)
    enddo

    print *, 'pblhmin ', pblhmin
    do k = 2, klev
       do ig = 1, klon
          if (teta(ig, 2).gt.teta(ig, 1)) then
             qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2
             kmin = kap*zlev(ig, k)*qmin
          else
             kmin = -1. ! kmin n'est utilise que pour les SL stables.
          endif
          if (kn(ig, k).lt.kmin.or.km(ig, k).lt.kmin) then
             kn(ig, k) = kmin
             km(ig, k) = kmin
             kq(ig, k) = kmin
             ! la longueur de melange est suposee etre l = kap z
             ! K = l q Sm d'ou q2 = (K/l Sm)**2
             q2(ig, k) = (qmin/sm(ig, k))**2
          endif
       enddo
    enddo

    ! Diagnostique pour stokage

    rino = rif
    smyam(:, 1:klev) = sm(:, 1:klev)
    styam = sm(:, 1:klev)*alpha(:, 1:klev)
    lyam(1:klon, 1:klev) = l(:, 1:klev)
    knyam(1:klon, 1:klev) = kn(:, 1:klev)

    first = .false.

  end SUBROUTINE yamada4

  !*******************************************************************

  function frif(ri)

    real frif
    real, intent(in):: ri

    frif = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))

  end function frif

  !*******************************************************************

  function falpha(ri)

    real falpha
    real, intent(in):: ri

    falpha = 1.318*(0.2231-ri)/(0.2341-ri)

  end function falpha

  !*******************************************************************

  function fsm(ri)

    real fsm
    real, intent(in):: ri

    fsm = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))

  end function fsm

  !*******************************************************************

  function fl(zzz, zl0, zq2, zn2)

    real fl
    real, intent(in):: zzz, zl0, zq2, zn2

    fl = max(min(zl0 * kap * zzz / (kap * zzz + zl0), &
         0.5 * sqrt(zq2) / sqrt(max(zn2, 1e-10))), 1.)

  end function fl

end module yamada4_m
1	guez	47	module yamada4_m
2	guez	3
3	guez	47	IMPLICIT NONE
4	guez	3
5	guez	47	real, parameter:: kap = 0.4
6			private
7			public yamada4
8	guez	3
9	guez	47	contains
10	guez	3
11	guez	47	SUBROUTINE yamada4(ngrid, dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, &
12			kq, ustar, iflag_pbl)
13	guez	3
14	guez	47	! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36
15	guez	3
16	guez	47	USE dimphy, ONLY : klev, klon
17	guez	3
18	guez	47	integer ngrid
19			REAL, intent(in):: dt ! pas de temps
20			real, intent(in):: g
21	guez	3
22	guez	47	REAL zlev(klon, klev+1)
23			! altitude à chaque niveau (interface inférieure de la couche de
24			! même indice)
25	guez	3
26	guez	47	REAL zlay(klon, klev) ! altitude au centre de chaque couche
27	guez	3
28	guez	47	REAL u(klon, klev), v(klon, klev)
29			! vitesse au centre de chaque couche (en entrée : la valeur au
30			! début du pas de temps)
31	guez	3
32	guez	47	REAL teta(klon, klev)
33			! température potentielle au centre de chaque couche (en entrée :
34			! la valeur au début du pas de temps)
35	guez	3
36	guez	47	REAL cd(klon) ! cdrag (en entrée : la valeur au début du pas de temps)
37	guez	3
38	guez	47	REAL, intent(inout):: q2(klon, klev+1)
39			! $q^2$ au bas de chaque couche
40			! En entrée : la valeur au début du pas de temps ; en sortie : la
41			! valeur à la fin du pas de temps.
42	guez	3
43	guez	47	REAL km(klon, klev+1)
44			! diffusivité turbulente de quantité de mouvement (au bas de
45			! chaque couche) (en sortie : la valeur à la fin du pas de temps)
46	guez	3
47	guez	47	REAL kn(klon, klev+1)
48			! diffusivité turbulente des scalaires (au bas de chaque couche)
49			! (en sortie : la valeur à la fin du pas de temps)
50	guez	3
51	guez	47	REAL kq(klon, klev+1)
52			real ustar(klon)
53	guez	3
54	guez	47	integer iflag_pbl
55			! iflag_pbl doit valoir entre 6 et 9
56			! l = 6, on prend systématiquement une longueur d'équilibre
57			! iflag_pbl = 6 : MY 2.0
58			! iflag_pbl = 7 : MY 2.0.Fournier
59			! iflag_pbl = 8 : MY 2.5
60			! iflag_pbl = 9 : un test ?
61	guez	3
62	guez	47	! Local:
63	guez	3
64	guez	47	real kmin, qmin, pblhmin(klon), coriol(klon)
65			real qpre
66			REAL unsdz(klon, klev)
67			REAL unsdzdec(klon, klev+1)
68			REAL kmpre(klon, klev+1), tmp2
69			REAL mpre(klon, klev+1)
70			real delta(klon, klev+1)
71			real aa(klon, klev+1), aa0, aa1
72			integer, PARAMETER:: nlay = klev
73			integer, PARAMETER:: nlev = klev+1
74			logical:: first = .true.
75			integer:: ipas = 0
76			integer ig, k
77			real ri
78			real rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
79			real m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
80			real dtetadz(klon, klev+1)
81			real m2cstat, mcstat, kmcstat
82			real l(klon, klev+1)
83			real, save:: l0(klon)
84			real sq(klon), sqz(klon), zz(klon, klev+1)
85			integer iter
86			real:: ric = 0.195, rifc = 0.191, b1 = 16.6, kap = 0.4
87			real rino(klon, klev+1), smyam(klon, klev), styam(klon, klev)
88			real lyam(klon, klev), knyam(klon, klev)
89	guez	3
90	guez	47	!-----------------------------------------------------------------------
91	guez	3
92	guez	47	if (.not. (iflag_pbl >= 6 .and. iflag_pbl <= 9)) then
93			print *, 'probleme de coherence dans appel a MY'
94			stop 1
95			endif
96	guez	3
97	guez	47	ipas = ipas+1
98	guez	3
99	guez	47	! les increments verticaux
100			DO ig = 1, ngrid
101			! alerte: zlev n'est pas declare a nlev
102			zlev(ig, nlev) = zlay(ig, nlay) +(zlay(ig, nlay) - zlev(ig, nlev-1))
103			ENDDO
104	guez	3
105	guez	47	DO k = 1, nlay
106			DO ig = 1, ngrid
107			unsdz(ig, k) = 1.E+0/(zlev(ig, k+1)-zlev(ig, k))
108			ENDDO
109			ENDDO
110			DO ig = 1, ngrid
111			unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))
112			ENDDO
113			DO k = 2, nlay
114			DO ig = 1, ngrid
115			unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))
116			ENDDO
117			ENDDO
118			DO ig = 1, ngrid
119			unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig, nlay+1)-zlay(ig, nlay))
120			ENDDO
121	guez	3
122	guez	47	do k = 2, klev
123			do ig = 1, ngrid
124			dz(ig, k) = zlay(ig, k)-zlay(ig, k-1)
125			m2(ig, k) = ((u(ig, k)-u(ig, k-1))2+(v(ig, k)-v(ig, k-1))2) &
126			/(dz(ig, k)*dz(ig, k))
127			dtetadz(ig, k) = (teta(ig, k)-teta(ig, k-1))/dz(ig, k)
128			n2(ig, k) = g2.dtetadz(ig, k)/(teta(ig, k-1)+teta(ig, k))
129			ri = n2(ig, k)/max(m2(ig, k), 1.e-10)
130			if (ri.lt.ric) then
131			rif(ig, k) = frif(ri)
132			else
133			rif(ig, k) = rifc
134			endif
135			if(rif(ig, k).lt.0.16) then
136			alpha(ig, k) = falpha(rif(ig, k))
137			sm(ig, k) = fsm(rif(ig, k))
138			else
139			alpha(ig, k) = 1.12
140			sm(ig, k) = 0.085
141			endif
142			zz(ig, k) = b1m2(ig, k)(1.-rif(ig, k))*sm(ig, k)
143			enddo
144			enddo
145	guez	3
146	guez	47	! Au premier appel, on détermine l et q2 de façon itérative.
147			! Itération pour déterminer la longueur de mélange
148	guez	3
149	guez	47	if (first .or. iflag_pbl == 6) then
150			do ig = 1, ngrid
151			l0(ig) = 10.
152			enddo
153			do k = 2, klev-1
154			do ig = 1, ngrid
155			l(ig, k) = l0(ig) * kap * zlev(ig, k) &
156			/ (kap * zlev(ig, k) + l0(ig))
157			enddo
158			enddo
159	guez	3
160	guez	47	do iter = 1, 10
161			do ig = 1, ngrid
162			sq(ig) = 1e-10
163			sqz(ig) = 1e-10
164			enddo
165			do k = 2, klev-1
166			do ig = 1, ngrid
167			q2(ig, k) = l(ig, k)*2 zz(ig, k)
168			l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
169			zq = sqrt(q2(ig, k))
170			sqz(ig) = sqz(ig) + zq * zlev(ig, k) &
171			* (zlay(ig, k) - zlay(ig, k-1))
172			sq(ig) = sq(ig) + zq * (zlay(ig, k) - zlay(ig, k-1))
173			enddo
174			enddo
175			do ig = 1, ngrid
176			l0(ig) = 0.2 * sqz(ig) / sq(ig)
177			enddo
178			enddo
179			endif
180	guez	3
181	guez	47	! Calcul de la longueur de melange.
182	guez	3
183	guez	47	! Mise a jour de l0
184			do ig = 1, ngrid
185			sq(ig) = 1.e-10
186			sqz(ig) = 1.e-10
187			enddo
188			do k = 2, klev-1
189			do ig = 1, ngrid
190			zq = sqrt(q2(ig, k))
191			sqz(ig) = sqz(ig)+zqzlev(ig, k)(zlay(ig, k)-zlay(ig, k-1))
192			sq(ig) = sq(ig)+zq*(zlay(ig, k)-zlay(ig, k-1))
193			enddo
194			enddo
195			do ig = 1, ngrid
196			l0(ig) = 0.2*sqz(ig)/sq(ig)
197			enddo
198			! calcul de l(z)
199			do k = 2, klev
200			do ig = 1, ngrid
201			l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
202			if(first) then
203			q2(ig, k) = l(ig, k)*2 zz(ig, k)
204			endif
205			enddo
206			enddo
207	guez	3
208	guez	47	! Yamada 2.0
209			if (iflag_pbl == 6) then
210			do k = 2, klev
211			do ig = 1, ngrid
212			q2(ig, k) = l(ig, k)*2 zz(ig, k)
213			enddo
214			enddo
215			else if (iflag_pbl == 7) then
216			! Yamada 2.Fournier
217	guez	3
218	guez	47	! Calcul de l, km, au pas precedent
219			do k = 2, klev
220			do ig = 1, ngrid
221			delta(ig, k) = q2(ig, k) / (l(ig, k)*2 sm(ig, k))
222			kmpre(ig, k) = l(ig, k) * sqrt(q2(ig, k)) * sm(ig, k)
223			mpre(ig, k) = sqrt(m2(ig, k))
224			enddo
225			enddo
226	guez	3
227	guez	47	do k = 2, klev-1
228			do ig = 1, ngrid
229			m2cstat = max(alpha(ig, k)*n2(ig, k)+delta(ig, k)/b1, 1.e-12)
230			mcstat = sqrt(m2cstat)
231	guez	3
232	guez	47	! puis on ecrit la valeur de q qui annule l'equation de m
233			! supposee en q3
234	guez	3
235	guez	47	IF (k == 2) THEN
236			kmcstat = 1.E+0 / mcstat &
237			(unsdz(ig, k)kmpre(ig, k+1) &
238			*mpre(ig, k+1) &
239			+unsdz(ig, k-1) &
240			*cd(ig) &
241			(sqrt(u(ig, 3)2+v(ig, 3)*2) &
242			-mcstat/unsdzdec(ig, k) &
243			-mpre(ig, k+1)/unsdzdec(ig, k+1))**2) &
244			/(unsdz(ig, k)+unsdz(ig, k-1))
245			ELSE
246			kmcstat = 1.E+0 / mcstat &
247			(unsdz(ig, k)kmpre(ig, k+1) &
248			*mpre(ig, k+1) &
249			+unsdz(ig, k-1)*kmpre(ig, k-1) &
250			*mpre(ig, k-1)) &
251			/(unsdz(ig, k)+unsdz(ig, k-1))
252			ENDIF
253			tmp2 = kmcstat / (sm(ig, k) / q2(ig, k)) /l(ig, k)
254			q2(ig, k) = max(tmp2, 1.e-12)**(2./3.)
255			enddo
256			enddo
257			else if (iflag_pbl >= 8) then
258			! Yamada 2.5 a la Didi
259	guez	3
260	guez	47	! Calcul de l, km, au pas precedent
261			do k = 2, klev
262			do ig = 1, ngrid
263			delta(ig, k) = q2(ig, k)/(l(ig, k)*2sm(ig, k))
264			if (delta(ig, k).lt.1.e-20) then
265			delta(ig, k) = 1.e-20
266			endif
267			km(ig, k) = l(ig, k)sqrt(q2(ig, k))sm(ig, k)
268			aa0 = (m2(ig, k)-alpha(ig, k)*n2(ig, k)-delta(ig, k)/b1)
269			aa1 = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1)
270			aa(ig, k) = aa1dt/(delta(ig, k)l(ig, k))
271			qpre = sqrt(q2(ig, k))
272			if (iflag_pbl == 8) then
273			if (aa(ig, k).gt.0.) then
274			q2(ig, k) = (qpre+aa(ig, k)qpreqpre)**2
275			else
276			q2(ig, k) = (qpre/(1.-aa(ig, k)qpre))*2
277			endif
278			else
279			! iflag_pbl = 9
280			if (aa(ig, k)*qpre.gt.0.9) then
281			q2(ig, k) = (qpre10.)*2
282			else
283			q2(ig, k) = (qpre/(1.-aa(ig, k)qpre))*2
284			endif
285			endif
286			q2(ig, k) = min(max(q2(ig, k), 1.e-10), 1.e4)
287			enddo
288			enddo
289			endif
290	guez	3
291	guez	47	! Calcul des coefficients de mélange
292			do k = 2, klev
293			do ig = 1, ngrid
294			zq = sqrt(q2(ig, k))
295			km(ig, k) = l(ig, k)zqsm(ig, k)
296			kn(ig, k) = km(ig, k)*alpha(ig, k)
297			kq(ig, k) = l(ig, k)zq0.2
298			enddo
299			enddo
300	guez	3
301	guez	47	! Traitement des cas noctrunes avec l'introduction d'une longueur
302			! minilale.
303	guez	3
304	guez	47	! Traitement particulier pour les cas tres stables.
305			! D'apres Holtslag Boville.
306	guez	3
307	guez	47	do ig = 1, ngrid
308			coriol(ig) = 1.e-4
309			pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546e-5)
310			enddo
311	guez	3
312	guez	47	print *, 'pblhmin ', pblhmin
313			do k = 2, klev
314			do ig = 1, klon
315			if (teta(ig, 2).gt.teta(ig, 1)) then
316			qmin = ustar(ig)(max(1.-zlev(ig, k)/pblhmin(ig), 0.))*2
317			kmin = kapzlev(ig, k)qmin
318			else
319			kmin = -1. ! kmin n'est utilise que pour les SL stables.
320			endif
321			if (kn(ig, k).lt.kmin.or.km(ig, k).lt.kmin) then
322			kn(ig, k) = kmin
323			km(ig, k) = kmin
324			kq(ig, k) = kmin
325			! la longueur de melange est suposee etre l = kap z
326			! K = l q Sm d'ou q2 = (K/l Sm)**2
327			q2(ig, k) = (qmin/sm(ig, k))**2
328			endif
329			enddo
330			enddo
331	guez	3
332	guez	47	! Diagnostique pour stokage
333	guez	3
334	guez	47	rino = rif
335			smyam(:, 1:klev) = sm(:, 1:klev)
336			styam = sm(:, 1:klev)*alpha(:, 1:klev)
337			lyam(1:klon, 1:klev) = l(:, 1:klev)
338			knyam(1:klon, 1:klev) = kn(:, 1:klev)
339	guez	3
340	guez	47	first = .false.
341	guez	3
342	guez	47	end SUBROUTINE yamada4
343	guez	3
344	guez	47	!*******************************************************************
345	guez	3
346	guez	47	function frif(ri)
347	guez	3
348	guez	47	real frif
349			real, intent(in):: ri
350	guez	3
351	guez	47	frif = 0.6588(ri+0.1776-sqrt(riri-0.3221*ri+0.03156))
352	guez	3
353	guez	47	end function frif
354	guez	3
355	guez	47	!*******************************************************************
356	guez	3
357	guez	47	function falpha(ri)
358	guez	3
359	guez	47	real falpha
360			real, intent(in):: ri
361	guez	3
362	guez	47	falpha = 1.318*(0.2231-ri)/(0.2341-ri)
363	guez	3
364	guez	47	end function falpha
365
366			!*******************************************************************
367
368			function fsm(ri)
369
370			real fsm
371			real, intent(in):: ri
372
373			fsm = 1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
374
375			end function fsm
376
377			!*******************************************************************
378
379			function fl(zzz, zl0, zq2, zn2)
380
381			real fl
382			real, intent(in):: zzz, zl0, zq2, zn2
383
384			fl = max(min(zl0 * kap * zzz / (kap * zzz + zl0), &
385			0.5 * sqrt(zq2) / sqrt(max(zn2, 1e-10))), 1.)
386
387			end function fl
388
389			end module yamada4_m