Changeset 2046
- Timestamp:
- May 18, 2014, 3:18:19 PM (11 years ago)
- Location:
- LMDZ5/trunk/libf/phylmd
- Files:
-
- 2 edited
-
1DUTILS.h (modified) (2 diffs)
-
thermcell_plume.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phylmd/1DUTILS.h
r2019 r2046 822 822 INTEGER np,ierr 823 823 REAL pi,x 824 REAL :: scaleheight=8. 824 825 825 826 !----------------------------------------------------------------------- … … 924 925 PRINT *, ap 925 926 926 DO l = 1, llm 927 dpres(l) = bp(l) - bp(l+1) 928 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 929 ENDDO 927 print*,'scaleheight=',scaleheight 928 DO l = 1, llm 929 dpres(l) = bp(l) - bp(l+1) 930 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 931 write(*, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & 932 log(preff/presnivs(l))*scaleheight & 933 , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ & 934 max(ap(l+1)+bp(l+1)*preff, 1.e-10)) 935 ENDDO 936 930 937 931 938 PRINT *,' PRESNIVS ' -
LMDZ5/trunk/libf/phylmd/thermcell_plume.F90
r2000 r2046 6 6 & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & 7 7 & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & 8 & ,lev_out,lunout1,igout) 8 & ,lev_out,lunout1,igout) 9 9 !-------------------------------------------------------------------------- 10 10 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance … … 80 80 real zbuoyjam(ngrid,klev) 81 81 real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis 82 real fact_shell 83 real ztv1,ztv2,factinv,zinv,zlmel 84 real ztv_est1,ztv_est2 82 85 real zcor,zdelta,zcvm5,qlbef 83 86 real betalpha,zbetalpha … … 92 95 ! Initialisation 93 96 94 ! print*,'THERMCELL PLUME OK'95 97 RLvCp = RLVTT/RCPD 96 98 fact_epsilon=0.002 97 99 betalpha=0.9 98 afact=2./3. 100 afact=2./3. 101 fact_shell=0.85 99 102 100 103 zbetalpha=betalpha/(1.+betalpha) … … 162 165 lalim(ig)=l+1 163 166 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) 164 ! print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)167 ! print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l) 165 168 endif 166 169 enddo … … 232 235 call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:)) 233 236 do ig=1,ngrid 234 ! print*,'active',active(ig),ig,l237 ! print*,'active',active(ig),ig,l 235 238 if(active(ig)) then 236 239 zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig)) … … 258 261 ! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2) 259 262 ! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact)) 263 264 !-------------------------------------------------- 265 !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l) 266 !-------------------------------------------------- 267 if (iflag_thermals_ed==8) then 268 ! Ancienne version 260 269 w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & 261 270 & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & 262 271 & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)) 263 if (w_est(ig,l+1).lt.0.) then 272 273 ! Nouvelle version Arnaud 274 else 275 w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & 276 & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & 277 & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) 278 endif 279 280 281 if (iflag_thermals_ed<6) then 282 zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) 283 ! fact_epsilon=0.0005/(zalpha+0.025)**0.5 284 ! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) 285 fact_epsilon=0.0002/(zalpha+0.1) 286 zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) 287 zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) 288 zdw2=afact*zbuoy(ig,l)/fact_epsilon 289 zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon 290 w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & 291 & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & 292 & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) 293 294 endif 295 !-------------------------------------------------- 296 !AJ052014: J'ai comment? ce if plus n?cessaire puisqu' 297 !on fait max(0.0001,.....) 298 !-------------------------------------------------- 299 300 ! if (w_est(ig,l+1).lt.0.) then 264 301 ! w_est(ig,l+1)=zw2(ig,l) 265 w_est(ig,l+1)=0.0001 266 endif 302 ! w_est(ig,l+1)=0.0001 303 ! endif 304 267 305 endif 268 306 enddo … … 296 334 297 335 lmel=fact_thermals_ed_dz*zlev(ig,l) 336 zlmel=zlev(ig,l)+lmel 298 337 ! lmel=2.5*(zlev(ig,l)-zlev(ig,l-1)) 299 338 lt=l+1 300 do it=1,klev-(l+1) 301 zdz2=zlev(ig,lt)-zlev(ig,l) 302 if (zdz2.gt.lmel) then 303 zdz3=zlev(ig,lt)-zlev(ig,lt-1) 339 zdz2=zlev(ig,lt)-zlev(ig,l) 340 !-------------------------------------------------- 341 !AJ052014: J'ai remplac? la boucle do par un do while 342 ! afin de faire moins de calcul dans la boucle 343 !-------------------------------------------------- 344 do while (zdz2.lt.lmel) 345 lt=lt+1 346 zdz2=zlev(ig,lt)-zlev(ig,l) 347 end do 348 349 zdz3=zlev(ig,lt)-zlev(ig,lt-1) 350 351 !-------------------------------------------------- 352 !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors 353 ! on cherche o? se trouve l'altitude d'inversion 354 ! en calculant ztv1 (interpolation de la valeur de 355 ! theta au niveau lt en utilisant les niveaux lt-1 et 356 ! lt-2) et ztv2 (interpolation avec les niveaux lt+1 357 ! et lt+2). Si theta r?ellement calcul?e au niveau lt 358 ! comprise entre ztv1 et ztv2, alors il y a inversion 359 ! et on calcule son altitude zinv en supposant que ztv(lt) 360 ! est une combinaison lin?aire de ztv1 et ztv2. 361 ! Ensuite, on calcule la flottabilit? en comparant 362 ! la temp?rature de la couche l ? celle de l'air situ? 363 ! l+lmel plus haut, ce qui n?cessite de savoir quel fraction 364 ! de cet air est au-dessus ou en-dessous de l'inversion 365 !-------------------------------------------------- 366 367 368 if (iflag_thermals_ed.lt.8) then 369 370 ztv1=(ztv(ig,lt-1)-ztv(ig,lt-2))*zlev(ig,lt)/(zlev(ig,lt-1)-zlev(ig,lt-2)) & 371 & +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) & 372 & /(zlev(ig,lt-1)-zlev(ig,lt-2)) 373 374 ztv2=(ztv(ig,lt+2)-ztv(ig,lt+1))*zlev(ig,lt)/(zlev(ig,lt+2)-zlev(ig,lt+1)) & 375 & +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) & 376 & /(zlev(ig,lt+2)-zlev(ig,lt+1)) 377 378 if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then 379 380 factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1) 381 zinv=zlev(ig,lt-1)+factinv*(zlev(ig,lt)-zlev(ig,lt-1)) 382 383 if (zlmel+0.5*zdz.ge.zinv) then 384 if (zlmel-0.5*zdz.ge.zinv) then 385 386 ztv_est(ig,l)=(ztv(ig,lt+2)-ztv(ig,lt+1))*(zlmel-0.*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) & 387 & +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) & 388 & /(zlev(ig,lt+2)-zlev(ig,lt+1)) 389 390 zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l) 391 392 else 393 394 ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))*0.5*(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) & 395 & +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) & 396 & /(zlev(ig,lt+2)-zlev(ig,lt+1)) 397 ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))*0.5*(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) & 398 & +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) & 399 & /(zlev(ig,lt-1)-zlev(ig,lt-2)) 400 zbuoyjam(ig,l)=fact_shell*RG*(((zlmel+0.5*zdz-zinv)/zdz)*(ztva_est(ig,l)- & 401 & ztv_est1)/ztv_est1+((zinv-zlmel+0.5*zdz)/zdz)*(ztva_est(ig,l)- & 402 & ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l) 403 404 endif 405 406 else 407 408 ztv_est(ig,l)=(ztv(ig,lt-1)-ztv(ig,lt-2))*(zlmel-0.*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) & 409 & +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) & 410 & /(zlev(ig,lt-1)-zlev(ig,lt-2)) 411 412 zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l) 413 ! ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))*0.5*(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) & 414 ! & +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) & 415 ! & /(zlev(ig,lt+2)-zlev(ig,lt+1)) 416 ! ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))*0.5*(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) & 417 ! & +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) & 418 ! & /(zlev(ig,lt-1)-zlev(ig,lt-2)) 419 ! zbuoyjam(ig,l)=fact_shell*RG*(((zlmel+0.5*zdz-zinv)/zdz)*(ztva_est(ig,l)- & 420 ! & ztv_est1)/ztv_est1+((zinv-zlmel+0.5*zdz)/zdz)*(ztva_est(ig,l)- & 421 ! & ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l) 422 423 424 425 426 ! print*,'on est pass? par l?',l,lt,zbuoyjam(ig,l),zbuoy(ig,l) 427 endif 428 429 430 else 431 304 432 ! ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l) 305 433 ! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) 434 435 zbuoyjam(ig,l)=fact_shell*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- & 436 & ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- & 437 & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l) 438 439 ! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- & 440 ! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- & 441 ! & po(ig,lt-1))/po(ig,lt-1)) 442 443 endif 444 445 else 306 446 307 447 zbuoyjam(ig,l)=1.*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- & … … 309 449 & ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l) 310 450 311 ! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- & 312 ! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- & 313 ! & po(ig,lt-1))/po(ig,lt-1)) 451 endif 314 452 315 else316 lt=lt+1317 endif318 enddo319 453 320 454 ! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) … … 325 459 ! entrbis=entr_star(ig,l) 326 460 327 328 detr_star(ig,l)=f_star(ig,l)*zdz & 461 if (iflag_thermals_ed.lt.6) then 462 fact_epsilon=0.0002/(zalpha+0.1) 463 endif 464 465 detr_star(ig,l)=f_star(ig,l)*zdz & 329 466 & *MAX(1.e-4, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m & 330 & + 0.012*(zdqt(ig,l)/zw2m)**0.5 467 & + 0.012*(zdqt(ig,l)/zw2m)**0.5) 331 468 332 469 zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) 333 470 334 471 entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & 335 & afact*zbuoy(ig,l)/zw2m - fact_epsilon ) 336 ! & afact*zbuoy(ig,l)/zw2m - fact_epsilon+ 0.012*(zdqt(ig,l)/zw2m)**0.5) 337 472 & afact*zbuoy(ig,l)/zw2m - fact_epsilon) 338 473 339 474 ! entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha* & … … 352 487 ! endif 353 488 354 ! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)489 ! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l) 355 490 ! Calcul du flux montant normalise 356 491 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & … … 395 530 zdzbis=zlev(ig,l+1)-zlev(ig,l-1) 396 531 zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) 397 532 fact_epsilon=0.002 398 533 zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) 399 534 zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) … … 404 539 & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & 405 540 & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) 541 542 if (iflag_thermals_ed.lt.6) then 543 zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l) 544 ! fact_epsilon=0.0005/(zalpha+0.025)**0.5 545 ! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) 546 fact_epsilon=0.0002/(zalpha+0.1)**1 547 zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) 548 zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) 549 zdw2= afact*zbuoy(ig,l)/(fact_epsilon) 550 zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon) 551 552 zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & 553 & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & 554 & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) 555 556 endif 557 558 406 559 endif 407 560 enddo … … 470 623 if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l 471 624 625 #undef wrgrads_thermcell 626 #ifdef wrgrads_thermcell 627 call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ') 628 call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ') 629 call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ') 630 call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ') 631 call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ') 632 call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ') 633 call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ') 634 #endif 635 636 472 637 return 473 638 end 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 474 670 475 671 … … 831 1027 zw2(ig,l+1)=0. 832 1028 linter(ig)=l+1 833 !CR:04/05/12:calcul linter 834 elseif (f_star(ig,l+1).lt.0.) then 835 linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & 836 & -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) 837 zw2(ig,l+1)=0. 838 !fin CR:04/05/12 839 endif 840 1029 endif 841 1030 842 1031 if (zw2(ig,l+1).lt.0.) then … … 881 1070 if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l 882 1071 1072 #undef wrgrads_thermcell 1073 #ifdef wrgrads_thermcell 1074 call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ') 1075 call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ') 1076 call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ') 1077 call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ') 1078 call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ') 1079 call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ') 1080 call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ') 1081 #endif 883 1082 884 1083
Note: See TracChangeset
for help on using the changeset viewer.
