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 à chaque niveau (interface inférieure de la couche de |
25 |
! même indice) |
! même 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ée : la valeur au |
31 |
! début du pas de temps) |
! début 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érature potentielle au centre de chaque couche (en entrée : |
35 |
! la valeur au début du pas de temps) |
! la valeur au début du pas de temps) |
36 |
|
|
37 |
REAL cd(klon) ! cdrag (en entrée : la valeur au début du pas de temps) |
REAL, intent(in):: cd(:) ! (ngrid) cdrag, valeur au début 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ée : la valeur au début du pas de temps ; en sortie : la |
42 |
! valeur à la fin du pas de temps. |
! valeur à 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é turbulente de quantité de mouvement (au bas de |
46 |
! chaque couche) (en sortie : la valeur à la fin du pas de temps) |
! chaque couche) (en sortie : la valeur à 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é turbulente des scalaires (au bas de chaque couche) |
50 |
! (en sortie : la valeur à la fin du pas de temps) |
! (en sortie : la valeur à 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ématiquement une longueur d'équilibre |
58 |
! iflag_pbl = 6 : MY 2.0 |
! iflag_pbl = 6 : MY 2.0 |
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), aa0, 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 |
87 |
real rino(klon, klev+1), smyam(klon, klev), styam(klon, klev) |
real rino(ngrid, klev+1), smyam(ngrid, klev), styam(ngrid, klev) |
88 |
real lyam(klon, klev), knyam(klon, klev) |
real lyam(ngrid, klev) |
89 |
|
|
90 |
!----------------------------------------------------------------------- |
!----------------------------------------------------------------------- |
91 |
|
|
92 |
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 |
|
93 |
|
|
94 |
ipas = ipas+1 |
ipas = ipas+1 |
95 |
|
|
96 |
! les increments verticaux |
! les increments verticaux |
97 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
98 |
! alerte: zlev n'est pas declare a nlev |
! alerte: zlev n'est pas declare a nlev |
99 |
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)) |
100 |
ENDDO |
ENDDO |
101 |
|
|
102 |
DO k = 1, nlay |
DO k = 1, klev |
103 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
104 |
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)) |
105 |
ENDDO |
ENDDO |
106 |
ENDDO |
ENDDO |
107 |
|
|
108 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
109 |
unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1)) |
unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1)) |
110 |
ENDDO |
ENDDO |
111 |
DO k = 2, nlay |
|
112 |
|
DO k = 2, klev |
113 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
114 |
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)) |
115 |
ENDDO |
ENDDO |
116 |
ENDDO |
ENDDO |
117 |
|
|
118 |
DO ig = 1, ngrid |
DO ig = 1, ngrid |
119 |
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)) |
120 |
ENDDO |
ENDDO |
121 |
|
|
122 |
do k = 2, klev |
do k = 2, klev |
132 |
else |
else |
133 |
rif(ig, k) = rifc |
rif(ig, k) = rifc |
134 |
endif |
endif |
135 |
if(rif(ig, k).lt.0.16) then |
if (rif(ig, k).lt.0.16) then |
136 |
alpha(ig, k) = falpha(rif(ig, k)) |
alpha(ig, k) = falpha(rif(ig, k)) |
137 |
sm(ig, k) = fsm(rif(ig, k)) |
sm(ig, k) = fsm(rif(ig, k)) |
138 |
else |
else |
199 |
do k = 2, klev |
do k = 2, klev |
200 |
do ig = 1, ngrid |
do ig = 1, ngrid |
201 |
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)) |
202 |
if(first) then |
if (first) then |
203 |
q2(ig, k) = l(ig, k)**2 * zz(ig, k) |
q2(ig, k) = l(ig, k)**2 * zz(ig, k) |
204 |
endif |
endif |
205 |
enddo |
enddo |
311 |
|
|
312 |
print *, 'pblhmin ', pblhmin |
print *, 'pblhmin ', pblhmin |
313 |
do k = 2, klev |
do k = 2, klev |
314 |
do ig = 1, klon |
do ig = 1, ngrid |
315 |
if (teta(ig, 2).gt.teta(ig, 1)) then |
if (teta(ig, 2).gt.teta(ig, 1)) then |
316 |
qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2 |
qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2 |
317 |
kmin = kap*zlev(ig, k)*qmin |
kmin = kap*zlev(ig, k)*qmin |
334 |
rino = rif |
rino = rif |
335 |
smyam(:, 1:klev) = sm(:, 1:klev) |
smyam(:, 1:klev) = sm(:, 1:klev) |
336 |
styam = sm(:, 1:klev)*alpha(:, 1:klev) |
styam = sm(:, 1:klev)*alpha(:, 1:klev) |
337 |
lyam(1:klon, 1:klev) = l(:, 1:klev) |
lyam(1:ngrid, 1:klev) = l(:, 1:klev) |
|
knyam(1:klon, 1:klev) = kn(:, 1:klev) |
|
338 |
|
|
339 |
first = .false. |
first = .false. |
340 |
|
|
342 |
|
|
343 |
!******************************************************************* |
!******************************************************************* |
344 |
|
|
345 |
function frif(ri) |
real function frif(ri) |
346 |
|
|
|
real frif |
|
347 |
real, intent(in):: ri |
real, intent(in):: ri |
348 |
|
|
349 |
frif = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) |
frif = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) |
352 |
|
|
353 |
!******************************************************************* |
!******************************************************************* |
354 |
|
|
355 |
function falpha(ri) |
real function falpha(ri) |
356 |
|
|
|
real falpha |
|
357 |
real, intent(in):: ri |
real, intent(in):: ri |
358 |
|
|
359 |
falpha = 1.318*(0.2231-ri)/(0.2341-ri) |
falpha = 1.318*(0.2231-ri)/(0.2341-ri) |
362 |
|
|
363 |
!******************************************************************* |
!******************************************************************* |
364 |
|
|
365 |
function fsm(ri) |
real function fsm(ri) |
366 |
|
|
|
real fsm |
|
367 |
real, intent(in):: ri |
real, intent(in):: ri |
368 |
|
|
369 |
fsm = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri)) |
fsm = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri)) |
372 |
|
|
373 |
!******************************************************************* |
!******************************************************************* |
374 |
|
|
375 |
function fl(zzz, zl0, zq2, zn2) |
real function fl(zzz, zl0, zq2, zn2) |
376 |
|
|
|
real fl |
|
377 |
real, intent(in):: zzz, zl0, zq2, zn2 |
real, intent(in):: zzz, zl0, zq2, zn2 |
378 |
|
|
379 |
fl = max(min(zl0 * kap * zzz / (kap * zzz + zl0), & |
fl = max(min(zl0 * kap * zzz / (kap * zzz + zl0), & |