2 |
|
|
3 |
IMPLICIT NONE |
IMPLICIT NONE |
4 |
|
|
|
real, parameter:: kap = 0.4 |
|
5 |
private |
private |
6 |
public yamada4 |
public yamada4 |
7 |
|
real, parameter:: kap = 0.4 |
8 |
|
|
9 |
contains |
contains |
10 |
|
|
11 |
SUBROUTINE yamada4(ngrid, dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, & |
SUBROUTINE yamada4(ngrid, dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, kq, & |
12 |
kq, ustar, iflag_pbl) |
ustar, iflag_pbl) |
13 |
|
|
14 |
! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36 |
! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36 |
15 |
|
|
16 |
USE dimphy, ONLY : klev, klon |
use nr_util, only: assert |
17 |
|
USE dimphy, ONLY: klev |
18 |
|
|
19 |
integer ngrid |
integer, intent(in):: ngrid |
20 |
REAL, intent(in):: dt ! pas de temps |
REAL, intent(in):: dt ! pas de temps |
21 |
real, intent(in):: g |
real, intent(in):: g |
22 |
|
|
23 |
REAL zlev(klon, klev+1) |
REAL zlev(ngrid, klev+1) |
24 |
! altitude à chaque niveau (interface inférieure de la couche de |
! altitude \`a chaque niveau (interface inf\'erieure de la couche de |
25 |
! même indice) |
! m\^eme indice) |
26 |
|
|
27 |
REAL zlay(klon, klev) ! altitude au centre de chaque couche |
REAL zlay(ngrid, klev) ! altitude au centre de chaque couche |
28 |
|
|
29 |
REAL u(klon, klev), v(klon, klev) |
REAL u(ngrid, klev), v(ngrid, klev) |
30 |
! vitesse au centre de chaque couche (en entrée : la valeur au |
! vitesse au centre de chaque couche (en entr\'ee : la valeur au |
31 |
! début du pas de temps) |
! d\'ebut du pas de temps) |
32 |
|
|
33 |
REAL teta(klon, klev) |
REAL, intent(in):: teta(ngrid, klev) |
34 |
! température potentielle au centre de chaque couche (en entrée : |
! temp\'erature potentielle au centre de chaque couche (en entr\'ee : |
35 |
! la valeur au début du pas de temps) |
! la valeur au d\'ebut du pas de temps) |
36 |
|
|
37 |
REAL, intent(in):: cd(:) ! (ngrid) cdrag, valeur au début du pas de temps |
REAL, intent(in):: cd(:) ! (ngrid) cdrag, valeur au d\'ebut du pas de temps |
38 |
|
|
39 |
REAL, intent(inout):: q2(klon, klev+1) |
REAL, intent(inout):: q2(ngrid, klev+1) |
40 |
! $q^2$ au bas de chaque couche |
! $q^2$ au bas de chaque couche |
41 |
! En entrée : la valeur au début du pas de temps ; en sortie : la |
! En entr\'ee : la valeur au d\'ebut du pas de temps ; en sortie : la |
42 |
! valeur à la fin du pas de temps. |
! valeur \`a la fin du pas de temps. |
43 |
|
|
44 |
REAL km(klon, klev+1) |
REAL km(ngrid, klev+1) |
45 |
! diffusivité turbulente de quantité de mouvement (au bas de |
! diffusivit\'e turbulente de quantit\'e de mouvement (au bas de |
46 |
! chaque couche) (en sortie : la valeur à la fin du pas de temps) |
! chaque couche) (en sortie : la valeur \`a la fin du pas de temps) |
47 |
|
|
48 |
REAL kn(klon, klev+1) |
REAL kn(ngrid, klev+1) |
49 |
! diffusivité turbulente des scalaires (au bas de chaque couche) |
! diffusivit\'e turbulente des scalaires (au bas de chaque couche) |
50 |
! (en sortie : la valeur à la fin du pas de temps) |
! (en sortie : la valeur \`a la fin du pas de temps) |
51 |
|
|
52 |
REAL kq(klon, klev+1) |
REAL kq(ngrid, klev+1) |
53 |
real ustar(klon) |
real ustar(ngrid) |
54 |
|
|
55 |
integer iflag_pbl |
integer, intent(in):: iflag_pbl |
56 |
! iflag_pbl doit valoir entre 6 et 9 |
! iflag_pbl doit valoir entre 6 et 9 |
57 |
! l = 6, on prend systématiquement une longueur d'équilibre |
! l = 6, on prend syst\'ematiquement une longueur d'\'equilibre |
58 |
! iflag_pbl = 6 : MY 2.0 |
! iflag_pbl = 6 : MY 2.0 |
59 |
! iflag_pbl = 7 : MY 2.0.Fournier |
! iflag_pbl = 7 : MY 2.0.Fournier |
60 |
! iflag_pbl = 8 : MY 2.5 |
! iflag_pbl = 8 : MY 2.5 |
62 |
|
|
63 |
! Local: |
! Local: |
64 |
|
|
65 |
real kmin, qmin, pblhmin(klon), coriol(klon) |
real kmin, qmin, pblhmin(ngrid), coriol(ngrid) |
66 |
real qpre |
real qpre |
67 |
REAL unsdz(klon, klev) |
REAL unsdz(ngrid, klev) |
68 |
REAL unsdzdec(klon, klev+1) |
REAL unsdzdec(ngrid, klev+1) |
69 |
REAL kmpre(klon, klev+1), tmp2 |
REAL kmpre(ngrid, klev+1), tmp2 |
70 |
REAL mpre(klon, klev+1) |
REAL mpre(ngrid, klev+1) |
71 |
real delta(klon, klev+1) |
real delta(ngrid, klev+1) |
72 |
real aa(klon, klev+1), aa0, aa1 |
real aa(ngrid, klev+1), aa1 |
|
integer, PARAMETER:: nlay = klev |
|
73 |
integer, PARAMETER:: nlev = klev+1 |
integer, PARAMETER:: nlev = klev+1 |
74 |
logical:: first = .true. |
logical:: first = .true. |
75 |
integer:: ipas = 0 |
integer:: ipas = 0 |
76 |
integer ig, k |
integer ig, k |
77 |
real ri |
real ri |
78 |
real rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev) |
real rif(ngrid, klev+1), sm(ngrid, klev+1), alpha(ngrid, klev) |
79 |
real m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1) |
real m2(ngrid, klev+1), dz(ngrid, klev+1), zq, n2(ngrid, klev+1) |
80 |
real dtetadz(klon, klev+1) |
real dtetadz(ngrid, klev+1) |
81 |
real m2cstat, mcstat, kmcstat |
real m2cstat, mcstat, kmcstat |
82 |
real l(klon, klev+1) |
real l(ngrid, klev+1) |
83 |
real, save:: l0(klon) |
real l0(ngrid) |
84 |
real sq(klon), sqz(klon), zz(klon, klev+1) |
real sq(ngrid), sqz(ngrid), zz(ngrid, klev+1) |
85 |
integer iter |
integer iter |
86 |
real:: ric = 0.195, rifc = 0.191, b1 = 16.6, kap = 0.4 |
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) |
|
87 |
|
|
88 |
!----------------------------------------------------------------------- |
!----------------------------------------------------------------------- |
89 |
|
|
90 |
if (.not. (iflag_pbl >= 6 .and. iflag_pbl <= 9)) then |
call assert(iflag_pbl >= 6 .and. iflag_pbl <= 9, "yamada4") |
|
print *, 'probleme de coherence dans appel a MY' |
|
|
stop 1 |
|
|
endif |
|
91 |
|
|
92 |
ipas = ipas+1 |
ipas = ipas+1 |
93 |
|
|
94 |
! les increments verticaux |
! les increments verticaux |
95 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
96 |
! alerte: zlev n'est pas declare a nlev |
! alerte: zlev n'est pas declare a nlev |
97 |
zlev(ig, nlev) = zlay(ig, nlay) +(zlay(ig, nlay) - zlev(ig, nlev-1)) |
zlev(ig, nlev) = zlay(ig, klev) +(zlay(ig, klev) - zlev(ig, nlev-1)) |
98 |
ENDDO |
ENDDO |
99 |
|
|
100 |
DO k = 1, nlay |
DO k = 1, klev |
101 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
102 |
unsdz(ig, k) = 1.E+0/(zlev(ig, k+1)-zlev(ig, k)) |
unsdz(ig, k) = 1.E+0/(zlev(ig, k+1)-zlev(ig, k)) |
103 |
ENDDO |
ENDDO |
104 |
ENDDO |
ENDDO |
105 |
|
|
106 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
107 |
unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1)) |
unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1)) |
108 |
ENDDO |
ENDDO |
109 |
DO k = 2, nlay |
|
110 |
|
DO k = 2, klev |
111 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
112 |
unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1)) |
unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1)) |
113 |
ENDDO |
ENDDO |
114 |
ENDDO |
ENDDO |
115 |
|
|
116 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
117 |
unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig, nlay+1)-zlay(ig, nlay)) |
unsdzdec(ig, klev+1) = 1.E+0/(zlev(ig, klev+1)-zlay(ig, klev)) |
118 |
ENDDO |
ENDDO |
119 |
|
|
120 |
do k = 2, klev |
do k = 2, klev |
130 |
else |
else |
131 |
rif(ig, k) = rifc |
rif(ig, k) = rifc |
132 |
endif |
endif |
133 |
if(rif(ig, k).lt.0.16) then |
if (rif(ig, k).lt.0.16) then |
134 |
alpha(ig, k) = falpha(rif(ig, k)) |
alpha(ig, k) = falpha(rif(ig, k)) |
135 |
sm(ig, k) = fsm(rif(ig, k)) |
sm(ig, k) = fsm(rif(ig, k)) |
136 |
else |
else |
141 |
enddo |
enddo |
142 |
enddo |
enddo |
143 |
|
|
144 |
! Au premier appel, on détermine l et q2 de façon itérative. |
! Au premier appel, on d\'etermine l et q2 de façon it\'erative. |
145 |
! Itération pour déterminer la longueur de mélange |
! It\'eration pour d\'eterminer la longueur de m\'elange |
146 |
|
|
147 |
if (first .or. iflag_pbl == 6) then |
if (first .or. iflag_pbl == 6) then |
148 |
do ig = 1, ngrid |
do ig = 1, ngrid |
197 |
do k = 2, klev |
do k = 2, klev |
198 |
do ig = 1, ngrid |
do ig = 1, ngrid |
199 |
l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k)) |
l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k)) |
200 |
if(first) then |
if (first) then |
201 |
q2(ig, k) = l(ig, k)**2 * zz(ig, k) |
q2(ig, k) = l(ig, k)**2 * zz(ig, k) |
202 |
endif |
endif |
203 |
enddo |
enddo |
263 |
delta(ig, k) = 1.e-20 |
delta(ig, k) = 1.e-20 |
264 |
endif |
endif |
265 |
km(ig, k) = l(ig, k)*sqrt(q2(ig, k))*sm(ig, k) |
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) |
|
266 |
aa1 = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1) |
aa1 = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1) |
267 |
aa(ig, k) = aa1*dt/(delta(ig, k)*l(ig, k)) |
aa(ig, k) = aa1*dt/(delta(ig, k)*l(ig, k)) |
268 |
qpre = sqrt(q2(ig, k)) |
qpre = sqrt(q2(ig, k)) |
285 |
enddo |
enddo |
286 |
endif |
endif |
287 |
|
|
288 |
! Calcul des coefficients de mélange |
! Calcul des coefficients de m\'elange |
289 |
do k = 2, klev |
do k = 2, klev |
290 |
do ig = 1, ngrid |
do ig = 1, ngrid |
291 |
zq = sqrt(q2(ig, k)) |
zq = sqrt(q2(ig, k)) |
308 |
|
|
309 |
print *, 'pblhmin ', pblhmin |
print *, 'pblhmin ', pblhmin |
310 |
do k = 2, klev |
do k = 2, klev |
311 |
do ig = 1, klon |
do ig = 1, ngrid |
312 |
if (teta(ig, 2).gt.teta(ig, 1)) then |
if (teta(ig, 2).gt.teta(ig, 1)) then |
313 |
qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2 |
qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2 |
314 |
kmin = kap*zlev(ig, k)*qmin |
kmin = kap*zlev(ig, k)*qmin |
326 |
enddo |
enddo |
327 |
enddo |
enddo |
328 |
|
|
|
! 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) |
|
|
|
|
329 |
first = .false. |
first = .false. |
330 |
|
|
331 |
end SUBROUTINE yamada4 |
end SUBROUTINE yamada4 |