phylmd/Interface_surf/yamada4.f

module yamada4_m

  IMPLICIT NONE

  private
  public yamada4
  real, parameter:: kap = 0.4

contains

  SUBROUTINE yamada4(zlev, zlay, u, v, teta, q2, coefm, coefh, ustar)

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

    ! Library:
    use nr_util, only: assert, assert_eq

    use comconst, only: dtphys
    USE conf_phys_m, ONLY: iflag_pbl
    USE dimphy, ONLY: klev
    USE suphec_m, ONLY: rg

    REAL, intent(in):: 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(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, intent(out):: coefm(:, 2:) ! (knon, 2:klev)
    ! 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, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
    ! diffusivit\'e turbulente des scalaires (au bas de chaque couche)
    ! (en sortie : la valeur \`a la fin du pas de temps)

    real, intent(in):: ustar(:) ! (knon)

    ! Local:
    
    integer knon
    real kmin, qmin
    real pblhmin(size(ustar)), coriol(size(ustar)) ! (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)
    logical:: first = .true.
    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(ustar)) ! (knon)
    real sq(size(ustar)), sqz(size(ustar)) ! (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(ustar), size(q2, 1), size(coefm, 1), &
         size(coefh, 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(coefm, 2) + 1, &
         size(coefh, 2) + 1], "yamada4 klev")

    ! les increments verticaux

    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) = rg*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, coefm, 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
             coefm(ig, k) = l(ig, k)*sqrt(q2(ig, k))*sm(ig, k)
             aa(ig, k) = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1)*dtphys/(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))
          coefm(ig, k) = l(ig, k)*zq*sm(ig, k)
          coefh(ig, k) = coefm(ig, k)*alpha(ig, k)
       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 (coefh(ig, k).lt.kmin.or.coefm(ig, k).lt.kmin) then
             coefh(ig, k) = kmin
             coefm(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

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

  pure 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

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

  pure real function falpha(ri)

    real, intent(in):: ri

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

  end function falpha

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

  pure 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

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

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