Sources/phylmd/yamada4.f

module yamada4_m

  IMPLICIT NONE

  private
  public yamada4
  real, parameter:: kap = 0.4

contains

  SUBROUTINE yamada4(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 nr_util, only: assert, assert_eq
    USE dimphy, ONLY: klev

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

    REAL zlev(:, :) ! (knon, klev + 1)
    ! altitude \`a chaque niveau (interface inf\'erieure de la couche de
    ! m\^eme indice)

    REAL, intent(in):: zlay(:, :) ! (knon, klev) altitude au centre de
                                  ! chaque couche

    REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev)
    ! vitesse au centre de chaque couche (en entr\'ee : la valeur au
    ! d\'ebut du pas de temps)

    REAL, intent(in):: teta(:, :) ! (knon, klev)
    ! temp\'erature potentielle au centre de chaque couche (en entr\'ee :
    ! la valeur au d\'ebut du pas de temps)

    REAL, intent(in):: cd(:) ! (knon) cdrag, valeur au d\'ebut du pas de temps

    REAL, intent(inout):: q2(:, :) ! (knon, klev + 1)
    ! $q^2$ au bas de chaque couche 
    ! En entr\'ee : la valeur au d\'ebut du pas de temps ; en sortie : la
    ! valeur \`a la fin du pas de temps.

    REAL km(:, :) ! (knon, klev + 1)
    ! diffusivit\'e turbulente de quantit\'e de mouvement (au bas de
    ! chaque couche) (en sortie : la valeur \`a la fin du pas de temps)

    REAL kn(:, :) ! (knon, klev + 1)
    ! diffusivit\'e turbulente des scalaires (au bas de chaque couche)
    ! (en sortie : la valeur \`a la fin du pas de temps)

    REAL kq(:, :) ! (knon, klev + 1)
    real, intent(in):: ustar(:) ! (knon)

    integer, intent(in):: iflag_pbl
    ! iflag_pbl doit valoir 6, 8 ou 9
    ! l = 6, on prend syst\'ematiquement une longueur d'\'equilibre
    ! iflag_pbl = 6 : Mellor and Yamada 2.0
    ! iflag_pbl = 8 : Mellor and Yamada 2.5
    ! iflag_pbl = 9 : un test ?

    ! Local:
    integer knon
    real kmin, qmin
    real pblhmin(size(cd)), coriol(size(cd)) ! (knon)
    real qpre
    REAL unsdz(size(zlay, 1), size(zlay, 2)) ! (knon, klev)
    REAL unsdzdec(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
    real delta(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
    real aa(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
    real aa1
    logical:: first = .true.
    integer:: ipas = 0
    integer ig, k
    real ri
    real, dimension(size(zlev, 1), size(zlev, 2)):: rif, sm ! (knon, klev + 1)
    real alpha(size(zlay, 1), size(zlay, 2)) ! (knon, klev)

    real, dimension(size(zlev, 1), size(zlev, 2)):: m2, dz, n2
    ! (knon, klev + 1)
    
    real zq
    real dtetadz(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
    real l(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
    real l0(size(cd)) ! (knon)
    real sq(size(cd)), sqz(size(cd)) ! (knon)
    real zz(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
    integer iter
    real:: ric = 0.195, rifc = 0.191, b1 = 16.6

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

    call assert(any(iflag_pbl == [6, 8, 9]), "yamada4 iflag_pbl")
    knon = assert_eq([size(zlev, 1), size(zlay, 1), size(u, 1), size(v, 1), &
         size(teta, 1), size(cd), size(q2, 1), size(km, 1), size(kn, 1), &
         size(kq, 1)], "yamada4 knon")
    call assert(klev == [size(zlev, 2) - 1, size(zlay, 2), size(u, 2), &
         size(v, 2), size(teta, 2), size(q2, 2) - 1, size(km, 2) - 1, &
         size(kn, 2) - 1, size(kq, 2) - 1], "yamada4 klev")

    ipas = ipas + 1

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

    DO k = 1, klev
       DO ig = 1, knon
          unsdz(ig, k) = 1.E+0/(zlev(ig, k + 1)-zlev(ig, k))
       ENDDO
    ENDDO

    DO ig = 1, knon
       unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))
    ENDDO

    DO k = 2, klev
       DO ig = 1, knon
          unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))
       ENDDO
    ENDDO

    DO ig = 1, knon
       unsdzdec(ig, klev + 1) = 1.E+0/(zlev(ig, klev + 1)-zlay(ig, klev))
    ENDDO

    do k = 2, klev
       do ig = 1, knon
          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\'etermine l et q2 de fa\ccon it\'erative.
    ! It\'eration pour d\'eterminer la longueur de m\'elange

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

       do iter = 1, 10
          do ig = 1, knon
             sq(ig) = 1e-10
             sqz(ig) = 1e-10
          enddo
          do k = 2, klev-1
             do ig = 1, knon
                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, knon
             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, knon
       sq(ig) = 1.e-10
       sqz(ig) = 1.e-10
    enddo
    do k = 2, klev-1
       do ig = 1, knon
          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, knon
       l0(ig) = 0.2*sqz(ig)/sq(ig)
    enddo
    ! calcul de l(z)
    do k = 2, klev
       do ig = 1, knon
          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

    if (iflag_pbl == 6) then
       ! Yamada 2.0
       do k = 2, klev
          do ig = 1, knon
             q2(ig, k) = l(ig, k)**2 * zz(ig, k)
          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, knon
             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)
             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\'elange
    do k = 2, klev
       do ig = 1, knon
          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, knon
       coriol(ig) = 1.e-4
       pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546e-5)
    enddo

    do k = 2, klev
       do ig = 1, knon
          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

    first = .false.

  end SUBROUTINE yamada4

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

  real function frif(ri)

    real, intent(in):: ri

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

  end function frif

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

  real function falpha(ri)

    real, intent(in):: ri

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

  end function falpha

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

  real function fsm(ri)

    real, intent(in):: ri

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

  end function fsm

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

  real function fl(zzz, zl0, zq2, zn2)

    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	118	private
6			public yamada4
7			real, parameter:: kap = 0.4
8	guez	3
9	guez	47	contains
10	guez	3
11	guez	227	SUBROUTINE yamada4(dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, kq, &
12	guez	118	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	227	use nr_util, only: assert, assert_eq
17	guez	118	USE dimphy, ONLY: klev
18	guez	3
19	guez	47	REAL, intent(in):: dt ! pas de temps
20			real, intent(in):: g
21	guez	3
22	guez	227	REAL zlev(:, :) ! (knon, klev + 1)
23	guez	178	! altitude \`a chaque niveau (interface inf\'erieure de la couche de
24			! m\^eme indice)
25	guez	3
26	guez	227	REAL, intent(in):: zlay(:, :) ! (knon, klev) altitude au centre de
27			! chaque couche
28	guez	3
29	guez	227	REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev)
30	guez	178	! vitesse au centre de chaque couche (en entr\'ee : la valeur au
31			! d\'ebut du pas de temps)
32	guez	3
33	guez	227	REAL, intent(in):: teta(:, :) ! (knon, klev)
34	guez	178	! temp\'erature potentielle au centre de chaque couche (en entr\'ee :
35			! la valeur au d\'ebut du pas de temps)
36	guez	3
37	guez	227	REAL, intent(in):: cd(:) ! (knon) cdrag, valeur au d\'ebut du pas de temps
38	guez	3
39	guez	227	REAL, intent(inout):: q2(:, :) ! (knon, klev + 1)
40	guez	47	! $q^2$ au bas de chaque couche
41	guez	178	! En entr\'ee : la valeur au d\'ebut du pas de temps ; en sortie : la
42			! valeur \`a la fin du pas de temps.
43	guez	3
44	guez	227	REAL km(:, :) ! (knon, klev + 1)
45	guez	178	! diffusivit\'e turbulente de quantit\'e de mouvement (au bas de
46			! chaque couche) (en sortie : la valeur \`a la fin du pas de temps)
47	guez	3
48	guez	227	REAL kn(:, :) ! (knon, klev + 1)
49	guez	178	! diffusivit\'e turbulente des scalaires (au bas de chaque couche)
50			! (en sortie : la valeur \`a la fin du pas de temps)
51	guez	3
52	guez	227	REAL kq(:, :) ! (knon, klev + 1)
53			real, intent(in):: ustar(:) ! (knon)
54	guez	3
55	guez	119	integer, intent(in):: iflag_pbl
56	guez	228	! iflag_pbl doit valoir 6, 8 ou 9
57	guez	178	! l = 6, on prend syst\'ematiquement une longueur d'\'equilibre
58	guez	228	! iflag_pbl = 6 : Mellor and Yamada 2.0
59			! iflag_pbl = 8 : Mellor and Yamada 2.5
60	guez	47	! iflag_pbl = 9 : un test ?
61	guez	3
62	guez	47	! Local:
63	guez	227	integer knon
64			real kmin, qmin
65			real pblhmin(size(cd)), coriol(size(cd)) ! (knon)
66	guez	47	real qpre
67	guez	227	REAL unsdz(size(zlay, 1), size(zlay, 2)) ! (knon, klev)
68			REAL unsdzdec(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
69			real delta(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
70			real aa(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
71			real aa1
72	guez	47	logical:: first = .true.
73			integer:: ipas = 0
74			integer ig, k
75			real ri
76	guez	227	real, dimension(size(zlev, 1), size(zlev, 2)):: rif, sm ! (knon, klev + 1)
77			real alpha(size(zlay, 1), size(zlay, 2)) ! (knon, klev)
78
79			real, dimension(size(zlev, 1), size(zlev, 2)):: m2, dz, n2
80			! (knon, klev + 1)
81
82			real zq
83			real dtetadz(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
84			real l(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
85			real l0(size(cd)) ! (knon)
86			real sq(size(cd)), sqz(size(cd)) ! (knon)
87			real zz(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
88	guez	47	integer iter
89	guez	196	real:: ric = 0.195, rifc = 0.191, b1 = 16.6
90	guez	3
91	guez	47	!-----------------------------------------------------------------------
92	guez	3
93	guez	228	call assert(any(iflag_pbl == [6, 8, 9]), "yamada4 iflag_pbl")
94	guez	227	knon = assert_eq([size(zlev, 1), size(zlay, 1), size(u, 1), size(v, 1), &
95			size(teta, 1), size(cd), size(q2, 1), size(km, 1), size(kn, 1), &
96			size(kq, 1)], "yamada4 knon")
97			call assert(klev == [size(zlev, 2) - 1, size(zlay, 2), size(u, 2), &
98			size(v, 2), size(teta, 2), size(q2, 2) - 1, size(km, 2) - 1, &
99			size(kn, 2) - 1, size(kq, 2) - 1], "yamada4 klev")
100	guez	3
101	guez	227	ipas = ipas + 1
102	guez	3
103	guez	47	! les increments verticaux
104	guez	227	DO ig = 1, knon
105			! alerte: zlev n'est pas declare a klev + 1
106			zlev(ig, klev + 1) = zlay(ig, klev) + (zlay(ig, klev) - zlev(ig, klev))
107	guez	47	ENDDO
108	guez	3
109	guez	118	DO k = 1, klev
110	guez	227	DO ig = 1, knon
111			unsdz(ig, k) = 1.E+0/(zlev(ig, k + 1)-zlev(ig, k))
112	guez	47	ENDDO
113			ENDDO
114	guez	118
115	guez	227	DO ig = 1, knon
116	guez	47	unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))
117			ENDDO
118	guez	118
119			DO k = 2, klev
120	guez	227	DO ig = 1, knon
121	guez	47	unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))
122			ENDDO
123			ENDDO
124	guez	118
125	guez	227	DO ig = 1, knon
126			unsdzdec(ig, klev + 1) = 1.E+0/(zlev(ig, klev + 1)-zlay(ig, klev))
127	guez	47	ENDDO
128	guez	3
129	guez	47	do k = 2, klev
130	guez	227	do ig = 1, knon
131	guez	47	dz(ig, k) = zlay(ig, k)-zlay(ig, k-1)
132	guez	227	m2(ig, k) = ((u(ig, k)-u(ig, k-1))2 + (v(ig, k)-v(ig, k-1))2) &
133	guez	47	/(dz(ig, k)*dz(ig, k))
134			dtetadz(ig, k) = (teta(ig, k)-teta(ig, k-1))/dz(ig, k)
135	guez	227	n2(ig, k) = g2.dtetadz(ig, k)/(teta(ig, k-1) + teta(ig, k))
136	guez	47	ri = n2(ig, k)/max(m2(ig, k), 1.e-10)
137			if (ri.lt.ric) then
138			rif(ig, k) = frif(ri)
139			else
140			rif(ig, k) = rifc
141			endif
142	guez	119	if (rif(ig, k).lt.0.16) then
143	guez	47	alpha(ig, k) = falpha(rif(ig, k))
144			sm(ig, k) = fsm(rif(ig, k))
145			else
146			alpha(ig, k) = 1.12
147			sm(ig, k) = 0.085
148			endif
149			zz(ig, k) = b1m2(ig, k)(1.-rif(ig, k))*sm(ig, k)
150			enddo
151			enddo
152	guez	3
153	guez	196	! Au premier appel, on d\'etermine l et q2 de fa\ccon it\'erative.
154	guez	178	! It\'eration pour d\'eterminer la longueur de m\'elange
155	guez	3
156	guez	47	if (first .or. iflag_pbl == 6) then
157	guez	227	do ig = 1, knon
158	guez	47	l0(ig) = 10.
159			enddo
160			do k = 2, klev-1
161	guez	227	do ig = 1, knon
162	guez	47	l(ig, k) = l0(ig) * kap * zlev(ig, k) &
163			/ (kap * zlev(ig, k) + l0(ig))
164			enddo
165			enddo
166	guez	3
167	guez	47	do iter = 1, 10
168	guez	227	do ig = 1, knon
169	guez	47	sq(ig) = 1e-10
170			sqz(ig) = 1e-10
171			enddo
172			do k = 2, klev-1
173	guez	227	do ig = 1, knon
174	guez	47	q2(ig, k) = l(ig, k)*2 zz(ig, k)
175			l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
176			zq = sqrt(q2(ig, k))
177			sqz(ig) = sqz(ig) + zq * zlev(ig, k) &
178			* (zlay(ig, k) - zlay(ig, k-1))
179			sq(ig) = sq(ig) + zq * (zlay(ig, k) - zlay(ig, k-1))
180			enddo
181			enddo
182	guez	227	do ig = 1, knon
183	guez	47	l0(ig) = 0.2 * sqz(ig) / sq(ig)
184			enddo
185			enddo
186			endif
187	guez	3
188	guez	47	! Calcul de la longueur de melange.
189	guez	3
190	guez	47	! Mise a jour de l0
191	guez	227	do ig = 1, knon
192	guez	47	sq(ig) = 1.e-10
193			sqz(ig) = 1.e-10
194			enddo
195			do k = 2, klev-1
196	guez	227	do ig = 1, knon
197	guez	47	zq = sqrt(q2(ig, k))
198	guez	227	sqz(ig) = sqz(ig) + zqzlev(ig, k)(zlay(ig, k)-zlay(ig, k-1))
199			sq(ig) = sq(ig) + zq*(zlay(ig, k)-zlay(ig, k-1))
200	guez	47	enddo
201			enddo
202	guez	227	do ig = 1, knon
203	guez	47	l0(ig) = 0.2*sqz(ig)/sq(ig)
204			enddo
205			! calcul de l(z)
206			do k = 2, klev
207	guez	227	do ig = 1, knon
208	guez	47	l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
209	guez	119	if (first) then
210	guez	47	q2(ig, k) = l(ig, k)*2 zz(ig, k)
211			endif
212			enddo
213			enddo
214	guez	3
215	guez	47	if (iflag_pbl == 6) then
216	guez	228	! Yamada 2.0
217	guez	47	do k = 2, klev
218	guez	227	do ig = 1, knon
219	guez	47	q2(ig, k) = l(ig, k)*2 zz(ig, k)
220			enddo
221			enddo
222			else if (iflag_pbl >= 8) then
223			! Yamada 2.5 a la Didi
224	guez	3
225	guez	47	! Calcul de l, km, au pas precedent
226			do k = 2, klev
227	guez	227	do ig = 1, knon
228	guez	47	delta(ig, k) = q2(ig, k)/(l(ig, k)*2sm(ig, k))
229			if (delta(ig, k).lt.1.e-20) then
230			delta(ig, k) = 1.e-20
231			endif
232			km(ig, k) = l(ig, k)sqrt(q2(ig, k))sm(ig, k)
233			aa1 = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1)
234			aa(ig, k) = aa1dt/(delta(ig, k)l(ig, k))
235			qpre = sqrt(q2(ig, k))
236			if (iflag_pbl == 8) then
237			if (aa(ig, k).gt.0.) then
238	guez	227	q2(ig, k) = (qpre + aa(ig, k)qpreqpre)**2
239	guez	47	else
240			q2(ig, k) = (qpre/(1.-aa(ig, k)qpre))*2
241			endif
242			else
243			! iflag_pbl = 9
244			if (aa(ig, k)*qpre.gt.0.9) then
245			q2(ig, k) = (qpre10.)*2
246			else
247			q2(ig, k) = (qpre/(1.-aa(ig, k)qpre))*2
248			endif
249			endif
250			q2(ig, k) = min(max(q2(ig, k), 1.e-10), 1.e4)
251			enddo
252			enddo
253			endif
254	guez	3
255	guez	178	! Calcul des coefficients de m\'elange
256	guez	47	do k = 2, klev
257	guez	227	do ig = 1, knon
258	guez	47	zq = sqrt(q2(ig, k))
259			km(ig, k) = l(ig, k)zqsm(ig, k)
260			kn(ig, k) = km(ig, k)*alpha(ig, k)
261			kq(ig, k) = l(ig, k)zq0.2
262			enddo
263			enddo
264	guez	3
265	guez	47	! Traitement des cas noctrunes avec l'introduction d'une longueur
266			! minilale.
267	guez	3
268	guez	47	! Traitement particulier pour les cas tres stables.
269			! D'apres Holtslag Boville.
270	guez	3
271	guez	227	do ig = 1, knon
272	guez	47	coriol(ig) = 1.e-4
273			pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546e-5)
274			enddo
275	guez	3
276	guez	47	do k = 2, klev
277	guez	227	do ig = 1, knon
278	guez	47	if (teta(ig, 2).gt.teta(ig, 1)) then
279			qmin = ustar(ig)(max(1.-zlev(ig, k)/pblhmin(ig), 0.))*2
280			kmin = kapzlev(ig, k)qmin
281			else
282			kmin = -1. ! kmin n'est utilise que pour les SL stables.
283			endif
284			if (kn(ig, k).lt.kmin.or.km(ig, k).lt.kmin) then
285			kn(ig, k) = kmin
286			km(ig, k) = kmin
287			kq(ig, k) = kmin
288			! la longueur de melange est suposee etre l = kap z
289			! K = l q Sm d'ou q2 = (K/l Sm)**2
290			q2(ig, k) = (qmin/sm(ig, k))**2
291			endif
292			enddo
293			enddo
294	guez	3
295	guez	47	first = .false.
296	guez	3
297	guez	47	end SUBROUTINE yamada4
298	guez	3
299	guez	47	!*******************************************************************
300	guez	3
301	guez	62	real function frif(ri)
302	guez	3
303	guez	47	real, intent(in):: ri
304	guez	3
305	guez	227	frif = 0.6588(ri + 0.1776-sqrt(riri-0.3221*ri + 0.03156))
306	guez	3
307	guez	47	end function frif
308	guez	3
309	guez	47	!*******************************************************************
310	guez	3
311	guez	62	real function falpha(ri)
312	guez	3
313	guez	47	real, intent(in):: ri
314	guez	3
315	guez	47	falpha = 1.318*(0.2231-ri)/(0.2341-ri)
316	guez	3
317	guez	47	end function falpha
318
319			!*******************************************************************
320
321	guez	62	real function fsm(ri)
322	guez	47
323			real, intent(in):: ri
324
325			fsm = 1.96(0.1912-ri)(0.2341-ri)/((1.-ri)*(0.2231-ri))
326
327			end function fsm
328
329			!*******************************************************************
330
331	guez	62	real function fl(zzz, zl0, zq2, zn2)
332	guez	47
333			real, intent(in):: zzz, zl0, zq2, zn2
334
335			fl = max(min(zl0 * kap * zzz / (kap * zzz + zl0), &
336			0.5 * sqrt(zq2) / sqrt(max(zn2, 1e-10))), 1.)
337
338			end function fl
339
340			end module yamada4_m