/[lmdze]/trunk/libf/phylmd/yamada4.f
ViewVC logotype

Contents of /trunk/libf/phylmd/yamada4.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (show annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 2 months ago) by guez
File size: 15671 byte(s)
Initial import
1 !
2 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/yamada4.F,v 1.1 2004/06/22 11:45:36 lmdzadmin Exp $
3 !
4 SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
5 s ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
6 s ,iflag_pbl)
7 use dimens_m
8 use dimphy
9 IMPLICIT NONE
10 c.......................................................................
11 c.......................................................................
12 c
13 c dt : pas de temps
14 c g : g
15 c zlev : altitude a chaque niveau (interface inferieure de la couche
16 c de meme indice)
17 c zlay : altitude au centre de chaque couche
18 c u,v : vitesse au centre de chaque couche
19 c (en entree : la valeur au debut du pas de temps)
20 c teta : temperature potentielle au centre de chaque couche
21 c (en entree : la valeur au debut du pas de temps)
22 c cd : cdrag
23 c (en entree : la valeur au debut du pas de temps)
24 c q2 : $q^2$ au bas de chaque couche
25 c (en entree : la valeur au debut du pas de temps)
26 c (en sortie : la valeur a la fin du pas de temps)
27 c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
28 c couche)
29 c (en sortie : la valeur a la fin du pas de temps)
30 c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
31 c (en sortie : la valeur a la fin du pas de temps)
32 c
33 c iflag_pbl doit valoir entre 6 et 9
34 c l=6, on prend systematiquement une longueur d'equilibre
35 c iflag_pbl=6 : MY 2.0
36 c iflag_pbl=7 : MY 2.0.Fournier
37 c iflag_pbl=8 : MY 2.5
38 c iflag_pbl=9 : un test ?
39
40 c.......................................................................
41 REAL dt,g,rconst
42 real plev(klon,klev+1),temp(klon,klev)
43 real ustar(klon)
44 real kmin,qmin,pblhmin(klon),coriol(klon)
45 REAL zlev(klon,klev+1)
46 REAL zlay(klon,klev)
47 REAL u(klon,klev)
48 REAL v(klon,klev)
49 REAL teta(klon,klev)
50 REAL cd(klon)
51 REAL q2(klon,klev+1),qpre
52 REAL unsdz(klon,klev)
53 REAL unsdzdec(klon,klev+1)
54
55 REAL km(klon,klev+1)
56 REAL kmpre(klon,klev+1),tmp2
57 REAL mpre(klon,klev+1)
58 REAL kn(klon,klev+1)
59 REAL kq(klon,klev+1)
60 real ff(klon,klev+1),delta(klon,klev+1)
61 real aa(klon,klev+1),aa0,aa1
62 integer iflag_pbl,ngrid
63
64
65 integer nlay,nlev
66 PARAMETER (nlay=klev)
67 PARAMETER (nlev=klev+1)
68
69 logical first
70 integer ipas
71 save first,ipas
72 data first,ipas/.true.,0/
73
74
75 integer ig,k
76
77
78 real ri,zrif,zalpha,zsm,zsn
79 real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
80
81 real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
82 real dtetadz(klon,klev+1)
83 real m2cstat,mcstat,kmcstat
84 real l(klon,klev+1),l0(klon)
85 save l0
86
87 real sq(klon),sqz(klon),zz(klon,klev+1)
88 integer iter
89
90 real ric,rifc,b1,kap
91 save ric,rifc,b1,kap
92 data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
93
94 real frif,falpha,fsm
95 real fl,zzz,zl0,zq2,zn2
96
97 real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
98 s ,lyam(klon,klev),knyam(klon,klev)
99 s ,w2yam(klon,klev),t2yam(klon,klev)
100 common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam
101
102 frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
103 falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
104 fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
105 fl(zzz,zl0,zq2,zn2)=
106 s max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
107 s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
108
109 if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
110 stop'probleme de coherence dans appel a MY'
111 endif
112
113 ipas=ipas+1
114 if (0.eq.1.and.first) then
115 do ig=1,1000
116 ri=(ig-800.)/500.
117 if (ri.lt.ric) then
118 zrif=frif(ri)
119 else
120 zrif=rifc
121 endif
122 if(zrif.lt.0.16) then
123 zalpha=falpha(zrif)
124 zsm=fsm(zrif)
125 else
126 zalpha=1.12
127 zsm=0.085
128 endif
129 c print*,ri,rif,zalpha,zsm
130 enddo
131 endif
132
133 c.......................................................................
134 c les increments verticaux
135 c.......................................................................
136 c
137 c!!!!! allerte !!!!!c
138 c!!!!! zlev n'est pas declare a nlev !!!!!c
139 c!!!!! ---->
140 DO ig=1,ngrid
141 zlev(ig,nlev)=zlay(ig,nlay)
142 & +( zlay(ig,nlay) - zlev(ig,nlev-1) )
143 ENDDO
144 c!!!!! <----
145 c!!!!! allerte !!!!!c
146 c
147 DO k=1,nlay
148 DO ig=1,ngrid
149 unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
150 ENDDO
151 ENDDO
152 DO ig=1,ngrid
153 unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
154 ENDDO
155 DO k=2,nlay
156 DO ig=1,ngrid
157 unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
158 ENDDO
159 ENDDO
160 DO ig=1,ngrid
161 unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
162 ENDDO
163 c
164 c.......................................................................
165
166 do k=2,klev
167 do ig=1,ngrid
168 dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
169 m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
170 s /(dz(ig,k)*dz(ig,k))
171 dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
172 n2(ig,k)=g*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
173 c n2(ig,k)=0.
174 ri=n2(ig,k)/max(m2(ig,k),1.e-10)
175 if (ri.lt.ric) then
176 rif(ig,k)=frif(ri)
177 else
178 rif(ig,k)=rifc
179 endif
180 if(rif(ig,k).lt.0.16) then
181 alpha(ig,k)=falpha(rif(ig,k))
182 sm(ig,k)=fsm(rif(ig,k))
183 else
184 alpha(ig,k)=1.12
185 sm(ig,k)=0.085
186 endif
187 zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
188 c print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
189
190
191 enddo
192 enddo
193
194
195 c====================================================================
196 c Au premier appel, on determine l et q2 de facon iterative.
197 c iterration pour determiner la longueur de melange
198
199
200 if (first.or.iflag_pbl.eq.6) then
201 do ig=1,ngrid
202 l0(ig)=10.
203 enddo
204 do k=2,klev-1
205 do ig=1,ngrid
206 l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
207 enddo
208 enddo
209
210 do iter=1,10
211 do ig=1,ngrid
212 sq(ig)=1.e-10
213 sqz(ig)=1.e-10
214 enddo
215 do k=2,klev-1
216 do ig=1,ngrid
217 q2(ig,k)=l(ig,k)**2*zz(ig,k)
218 l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
219 zq=sqrt(q2(ig,k))
220 sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
221 sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
222 enddo
223 enddo
224 do ig=1,ngrid
225 l0(ig)=0.2*sqz(ig)/sq(ig)
226 c l0(ig)=30.
227 enddo
228 c print*,'ITER=',iter,' L0=',l0
229
230 enddo
231
232 c print*,'Fin de l initialisation de q2 et l0'
233
234 endif ! first
235
236 c====================================================================
237 c Calcul de la longueur de melange.
238 c====================================================================
239
240 c Mise a jour de l0
241 do ig=1,ngrid
242 sq(ig)=1.e-10
243 sqz(ig)=1.e-10
244 enddo
245 do k=2,klev-1
246 do ig=1,ngrid
247 zq=sqrt(q2(ig,k))
248 sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
249 sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
250 enddo
251 enddo
252 do ig=1,ngrid
253 l0(ig)=0.2*sqz(ig)/sq(ig)
254 c l0(ig)=30.
255 enddo
256 c print*,'ITER=',iter,' L0=',l0
257 c calcul de l(z)
258 do k=2,klev
259 do ig=1,ngrid
260 l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
261 if(first) then
262 q2(ig,k)=l(ig,k)**2*zz(ig,k)
263 endif
264 enddo
265 enddo
266
267 c====================================================================
268 c Yamada 2.0
269 c====================================================================
270 if (iflag_pbl.eq.6) then
271
272 do k=2,klev
273 do ig=1,ngrid
274 q2(ig,k)=l(ig,k)**2*zz(ig,k)
275 enddo
276 enddo
277
278
279 else if (iflag_pbl.eq.7) then
280 c====================================================================
281 c Yamada 2.Fournier
282 c====================================================================
283
284 c Calcul de l, km, au pas precedent
285 do k=2,klev
286 do ig=1,ngrid
287 c print*,'SMML=',sm(ig,k),l(ig,k)
288 delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
289 kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
290 mpre(ig,k)=sqrt(m2(ig,k))
291 c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
292 enddo
293 enddo
294
295 do k=2,klev-1
296 do ig=1,ngrid
297 m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
298 mcstat=sqrt(m2cstat)
299
300 c print*,'M2 L=',k,mpre(ig,k),mcstat
301 c
302 c -----{puis on ecrit la valeur de q qui annule l'equation de m
303 c supposee en q3}
304 c
305 IF (k.eq.2) THEN
306 kmcstat=1.E+0 / mcstat
307 & *( unsdz(ig,k)*kmpre(ig,k+1)
308 & *mpre(ig,k+1)
309 & +unsdz(ig,k-1)
310 & *cd(ig)
311 & *( sqrt(u(ig,3)**2+v(ig,3)**2)
312 & -mcstat/unsdzdec(ig,k)
313 & -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
314 & /( unsdz(ig,k)+unsdz(ig,k-1) )
315 ELSE
316 kmcstat=1.E+0 / mcstat
317 & *( unsdz(ig,k)*kmpre(ig,k+1)
318 & *mpre(ig,k+1)
319 & +unsdz(ig,k-1)*kmpre(ig,k-1)
320 & *mpre(ig,k-1) )
321 & /( unsdz(ig,k)+unsdz(ig,k-1) )
322 ENDIF
323 c print*,'T2 L=',k,tmp2
324 tmp2=kmcstat
325 & /( sm(ig,k)/q2(ig,k) )
326 & /l(ig,k)
327 q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
328 c print*,'Q2 L=',k,q2(ig,k)
329 c
330 enddo
331 enddo
332
333 else if (iflag_pbl.ge.8) then
334 c====================================================================
335 c Yamada 2.5 a la Didi
336 c====================================================================
337
338
339 c Calcul de l, km, au pas precedent
340 do k=2,klev
341 do ig=1,ngrid
342 c print*,'SMML=',sm(ig,k),l(ig,k)
343 delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
344 if (delta(ig,k).lt.1.e-20) then
345 c print*,'ATTENTION L=',k,' Delta=',delta(ig,k)
346 delta(ig,k)=1.e-20
347 endif
348 km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
349 aa0=
350 s (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
351 aa1=
352 s (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
353 c abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
354 aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))
355 c print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
356 qpre=sqrt(q2(ig,k))
357 if (iflag_pbl.eq.8 ) then
358 if (aa(ig,k).gt.0.) then
359 q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
360 else
361 q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
362 endif
363 else ! iflag_pbl=9
364 if (aa(ig,k)*qpre.gt.0.9) then
365 q2(ig,k)=(qpre*10.)**2
366 else
367 q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
368 endif
369 endif
370 q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
371 c print*,'Q2 L=',k,q2(ig,k),qpre*qpre
372 enddo
373 enddo
374
375 endif ! Fin du cas 8
376
377 c print*,'OK8'
378
379 c====================================================================
380 c Calcul des coefficients de mélange
381 c====================================================================
382 do k=2,klev
383 c print*,'k=',k
384 do ig=1,ngrid
385 cabde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
386 zq=sqrt(q2(ig,k))
387 km(ig,k)=l(ig,k)*zq*sm(ig,k)
388 kn(ig,k)=km(ig,k)*alpha(ig,k)
389 kq(ig,k)=l(ig,k)*zq*0.2
390 c print*,'KML=',km(ig,k),kn(ig,k)
391 enddo
392 enddo
393
394 c if (iflag_pbl.ge.7..and.0.eq.1) then
395 c q2(:,1)=q2(:,2)
396 c call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
397 c endif
398
399 c Traitement des cas noctrunes avec l'introduction d'une longueur
400 c minilale.
401
402 c====================================================================
403 c Traitement particulier pour les cas tres stables.
404 c D'apres Holtslag Boville.
405
406 print*,'YAMADA4 0'
407
408 do ig=1,ngrid
409 coriol(ig)=1.e-4
410 pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
411 enddo
412
413 print*,'pblhmin ',pblhmin
414 CTest a remettre 21 11 02
415 c test abd 13 05 02 if(0.eq.1) then
416 if(1.eq.1) then
417 do k=2,klev
418 do ig=1,klon
419 if (teta(ig,2).gt.teta(ig,1)) then
420 qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
421 kmin=kap*zlev(ig,k)*qmin
422 else
423 kmin=-1. ! kmin n'est utilise que pour les SL stables.
424 endif
425 if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
426 c print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
427 c s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
428 kn(ig,k)=kmin
429 km(ig,k)=kmin
430 kq(ig,k)=kmin
431 c la longueur de melange est suposee etre l= kap z
432 c K=l q Sm d'ou q2=(K/l Sm)**2
433 q2(ig,k)=(qmin/sm(ig,k))**2
434 endif
435 enddo
436 enddo
437 endif
438
439 print*,'YAMADA4 1'
440 c Diagnostique pour stokage
441
442 rino=rif
443 smyam(:,1:klev)=sm(:,1:klev)
444 styam=sm(:,1:klev)*alpha(:,1:klev)
445 lyam(1:klon,1:klev)=l(:,1:klev)
446 knyam(1:klon,1:klev)=kn(:,1:klev)
447
448 c Estimations de w'2 et T'2 d'apres Abdela et McFarlane
449
450 if(1.eq.0)then
451 w2yam=q2(:,1:klev)*0.24
452 s +lyam(:,1:klev)*5.17*kn(:,1:klev)*n2(:,1:klev)
453 s /sqrt(q2(:,1:klev))
454
455 t2yam=9.1*kn(:,1:klev)*dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))
456 s *lyam(:,1:klev)
457 endif
458
459 c print*,'OKFIN'
460 first=.false.
461 return
462 end

  ViewVC Help
Powered by ViewVC 1.1.21