Sources/phylmd/yamada4.f

module yamada4_m

  IMPLICIT NONE

  private
  public yamada4
  real, parameter:: kap = 0.4

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

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

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

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

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

    REAL, intent(in):: teta(ngrid, 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(:) ! (ngrid) cdrag, valeur au d\'ebut du pas de temps

    REAL, intent(inout):: q2(ngrid, 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(ngrid, 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(ngrid, 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(ngrid, klev+1)
    real ustar(ngrid)

    integer, intent(in):: iflag_pbl
    ! iflag_pbl doit valoir entre 6 et 9
    ! l = 6, on prend syst\'ematiquement une longueur d'\'equilibre
    ! 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(ngrid), coriol(ngrid)
    real qpre
    REAL unsdz(ngrid, klev)
    REAL unsdzdec(ngrid, klev+1)
    REAL kmpre(ngrid, klev+1), tmp2
    REAL mpre(ngrid, klev+1)
    real delta(ngrid, klev+1)
    real aa(ngrid, klev+1), aa1
    integer, PARAMETER:: nlev = klev+1
    logical:: first = .true.
    integer:: ipas = 0
    integer ig, k
    real ri
    real rif(ngrid, klev+1), sm(ngrid, klev+1), alpha(ngrid, klev)
    real m2(ngrid, klev+1), dz(ngrid, klev+1), zq, n2(ngrid, klev+1)
    real dtetadz(ngrid, klev+1)
    real m2cstat, mcstat, kmcstat
    real l(ngrid, klev+1)
    real l0(ngrid)
    real sq(ngrid), sqz(ngrid), zz(ngrid, klev+1)
    integer iter
    real:: ric = 0.195, rifc = 0.191, b1 = 16.6

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

    call assert(iflag_pbl >= 6 .and. iflag_pbl <= 9, "yamada4")

    ipas = ipas+1

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

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