- Timestamp:
- Feb 18, 2010, 2:14:02 PM (14 years ago)
- Location:
- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd
- Files:
-
- 10 edited
-
calltherm.F90 (modified) (2 diffs)
-
pbl_surface_mod.F90 (modified) (1 diff)
-
phyetat0.F (modified) (2 diffs)
-
physiq.F (modified) (2 diffs)
-
printflag.F (modified) (1 diff)
-
thermcell_closure.F90 (modified) (1 diff)
-
thermcell_height.F90 (modified) (1 diff)
-
thermcell_main.F90 (modified) (6 diffs)
-
thermcell_plume.F90 (modified) (2 diffs)
-
yamada4.F (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90
r1299 r1311 204 204 & ,zmax0,f0,zw2,fraca) 205 205 else if (iflag_thermals==15.or.iflag_thermals==16) then 206 207 ! print*,'THERM iflag_thermas_ed=',iflag_thermals_ed 206 208 CALL thermcell_main(itap,klon,klev,zdt & 207 209 & ,pplay,paprs,pphi,debut & … … 210 212 & ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax & 211 213 & ,ratqscth,ratqsdiff,zqsatth & 212 & ,r_aspect_thermals,l_mix_thermals 213 & ,tau_thermals, Ale,Alp,lalim_conv,wght_th &214 & ,r_aspect_thermals,l_mix_thermals & 215 & ,tau_thermals,iflag_thermals_ed,Ale,Alp,lalim_conv,wght_th & 214 216 & ,zmax0,f0,zw2,fraca) 215 217 if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK' -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/pbl_surface_mod.F90
r1299 r1311 484 484 485 485 ! Initialize ok_flux_surf (for 1D model) 486 ok_flux_surf=.FALSE.486 if (klon>1) ok_flux_surf=.FALSE. 487 487 488 488 ! Initilize debug IO -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyetat0.F
r1298 r1311 228 228 $ 'coherente ', i, zmasq(i), pctsrf(i, is_ter) 229 229 $ ,pctsrf(i, is_lic) 230 zmasq(i) = fractint(i) 230 231 ENDIF 231 232 END DO … … 237 238 $ 'coherente ', i, zmasq(i) , pctsrf(i, is_oce) 238 239 $ ,pctsrf(i, is_sic) 240 zmasq(i) = fractint(i) 239 241 ENDIF 240 242 END DO -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F
r1299 r1311 2339 2339 con rajoute ale et alp, et les caracteristiques de la couche alim 2340 2340 s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca) 2341 2342 ! ---------------------------------------------------------------------- 2343 ! Transport de la TKE par les panaches thermiques. 2344 ! FH : 2010/02/01 2345 if (iflag_pbl.eq.10) then 2346 call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm, 2347 s rg,paprs,pbl_tke) 2348 endif 2349 ! ---------------------------------------------------------------------- 2350 2341 2351 endif 2352 2342 2353 2343 2354 … … 2550 2561 2551 2562 ! les ratqs sont une combinaison de ratqss et ratqsc 2552 write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs2563 ! write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs 2553 2564 2554 2565 if (tau_ratqs>1.e-10) then -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/printflag.F
r1279 r1311 85 85 IF( INT( tabcntr0( 6 ) ) .NE. nbapp_rad ) THEN 86 86 PRINT 21, INT(tabcntr0(6)), nbapp_rad 87 radpas0 = NINT( 86400./tabcntr0(1)/INT( tabcntr0(6) ) )87 ! radpas0 = NINT( 86400./tabcntr0(1)/INT( tabcntr0(6) ) ) 88 88 PRINT 100 89 89 PRINT 22, radpas0, radpas -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_closure.F90
r1294 r1311 39 39 40 40 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 41 print*,'THERMCELL CLOSURE 26E'41 !print*,'THERMCELL CLOSURE 26E' 42 42 43 43 alim_star2(:)=0. -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_height.F90
r1026 r1311 40 40 enddo 41 41 enddo 42 43 ! On traite le cas particulier qu'il faudrait éviter ou le thermique 44 ! atteind le haut du modele ... 45 do ig=1,ngrid 46 if ( zw2(ig,nlay) > 1.e-10 ) then 47 print*,'WARNING !!!!! W2 thermiques non nul derniere couche ' 48 lmax(ig)=nlay 49 endif 50 enddo 51 42 52 ! pas de thermique si couche 1 stable 43 53 do ig=1,ngrid -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90
r1299 r1311 8 8 & ,fm0,entr0,detr0,zqta,zqla,lmax & 9 9 & ,ratqscth,ratqsdiff,zqsatth & 10 & ,r_aspect,l_mix,tau_thermals &10 & ,r_aspect,l_mix,tau_thermals,iflag_thermals_ed & 11 11 & ,Ale_bl,Alp_bl,lalim_conv,wght_th & 12 12 & ,zmax0, f0,zw2,fraca) … … 55 55 INTEGER ngrid,nlay,w2di 56 56 real tau_thermals 57 integer iflag_thermals_ed 57 58 real ptimestep,l_mix,r_aspect 58 59 REAL pt(ngrid,nlay),pdtadj(ngrid,nlay) … … 228 229 ENDIF 229 230 ! 230 write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'231 ! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)' 231 232 do ig=1,klon 232 233 if (prt_level.ge.20) then … … 380 381 if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out 381 382 !IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & 382 CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & 383 & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & 384 & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & 385 & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & 386 & ,lev_out,lunout1,igout) 383 384 ! Gestion temporaire de plusieurs appels à thermcell_plume au travers 385 ! de la variable iflag_thermals 386 387 ! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed 388 if (iflag_thermals_ed<=9) then 389 ! print*,'THERM NOVUELLE/NOUVELLE/ANCIENNE' 390 CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,& 391 & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & 392 & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & 393 & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & 394 & ,lev_out,lunout1,igout) 395 396 elseif (iflag_thermals_ed<=19) then 397 ! Version d'Arnaud Jam 398 ! print*,'THERM RIO et al 2010' 399 CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,& 400 & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & 401 & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & 402 & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & 403 & ,lev_out,lunout1,igout) 404 405 endif 406 387 407 if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out 388 408 … … 421 441 ! 422 442 ! 423 write(lunout,*)'THERM NOUVEAU RIO2009 31B'443 !! write(lunout,*)'THERM NOUVEAU XXXXX' 424 444 CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & 425 445 & lalim,lmin,zmax_sec,wmax_sec,lev_out) … … 440 460 441 461 442 print*,'THERM 26JJJ'462 !print*,'THERM 26JJJ' 443 463 444 464 ! Choix de la fonction d'alimentation utilisee pour la fermeture. -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90
r1299 r1311 114 114 endif 115 115 116 write(lunout,*)'THERM 31H '116 ! write(lunout,*)'THERM 31H ' 117 117 118 118 ! Initialisations des variables reeles … … 483 483 return 484 484 end 485 486 487 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 488 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 489 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 490 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 491 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 492 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 493 SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, & 494 & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & 495 & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & 496 & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & 497 & ,lev_out,lunout1,igout) 498 499 !-------------------------------------------------------------------------- 500 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance 501 ! Version conforme a l'article de Rio et al. 2010. 502 ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin 503 !-------------------------------------------------------------------------- 504 505 IMPLICIT NONE 506 507 #include "YOMCST.h" 508 #include "YOETHF.h" 509 #include "FCTTRE.h" 510 #include "iniprint.h" 511 #include "thermcell.h" 512 513 INTEGER itap 514 INTEGER lunout1,igout 515 INTEGER ngrid,klev 516 REAL ptimestep 517 REAL ztv(ngrid,klev) 518 REAL zthl(ngrid,klev) 519 REAL po(ngrid,klev) 520 REAL zl(ngrid,klev) 521 REAL rhobarz(ngrid,klev) 522 REAL zlev(ngrid,klev+1) 523 REAL pplev(ngrid,klev+1) 524 REAL pphi(ngrid,klev) 525 REAL zpspsk(ngrid,klev) 526 REAL alim_star(ngrid,klev) 527 REAL f0(ngrid) 528 INTEGER lalim(ngrid) 529 integer lev_out ! niveau pour les print 530 531 real alim_star_tot(ngrid) 532 533 REAL ztva(ngrid,klev) 534 REAL ztla(ngrid,klev) 535 REAL zqla(ngrid,klev) 536 REAL zqta(ngrid,klev) 537 REAL zha(ngrid,klev) 538 539 REAL detr_star(ngrid,klev) 540 REAL coefc 541 REAL entr_star(ngrid,klev) 542 REAL detr(ngrid,klev) 543 REAL entr(ngrid,klev) 544 545 REAL csc(ngrid,klev) 546 547 REAL zw2(ngrid,klev+1) 548 REAL w_est(ngrid,klev+1) 549 REAL f_star(ngrid,klev+1) 550 REAL wa_moy(ngrid,klev+1) 551 552 REAL ztva_est(ngrid,klev) 553 REAL zqla_est(ngrid,klev) 554 REAL zqsatth(ngrid,klev) 555 REAL zta_est(ngrid,klev) 556 REAL zdw2 557 REAL zw2modif 558 REAL zw2fact 559 REAL zeps(ngrid,klev) 560 561 REAL linter(ngrid) 562 INTEGER lmix(ngrid) 563 INTEGER lmix_bis(ngrid) 564 REAL wmaxa(ngrid) 565 566 INTEGER ig,l,k 567 568 real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m 569 real zbuoybis 570 real zcor,zdelta,zcvm5,qlbef,zdz2 571 real betalpha,zbetalpha 572 real Tbef,qsatbef,b1,eps, afact 573 real dqsat_dT,DT,num,denom,m 574 REAL REPS,RLvCp,DDT0 575 PARAMETER (DDT0=.01) 576 logical Zsat 577 LOGICAL active(ngrid),activetmp(ngrid) 578 REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2 579 REAL c2(ngrid,klev) 580 Zsat=.false. 581 ! Initialisation 582 583 RLvCp = RLVTT/RCPD 584 fact_epsilon=0.002 585 betalpha=0.9 586 afact=2./3. 587 588 zbetalpha=betalpha/(1.+betalpha) 589 590 ! print*,'THERM 31B' 591 592 ! Initialisations des variables reeles 593 if (1==0) then 594 ztva(:,:)=ztv(:,:) 595 ztva_est(:,:)=ztva(:,:) 596 ztla(:,:)=zthl(:,:) 597 zqta(:,:)=po(:,:) 598 zha(:,:) = ztva(:,:) 599 else 600 ztva(:,:)=0. 601 ztva_est(:,:)=0. 602 ztla(:,:)=0. 603 zqta(:,:)=0. 604 zha(:,:) =0. 605 endif 606 607 zqla_est(:,:)=0. 608 zqsatth(:,:)=0. 609 zqla(:,:)=0. 610 detr_star(:,:)=0. 611 entr_star(:,:)=0. 612 alim_star(:,:)=0. 613 alim_star_tot(:)=0. 614 csc(:,:)=0. 615 detr(:,:)=0. 616 entr(:,:)=0. 617 zw2(:,:)=0. 618 zbuoy(:,:)=0. 619 gamma(:,:)=0. 620 zeps(:,:)=0. 621 w_est(:,:)=0. 622 f_star(:,:)=0. 623 wa_moy(:,:)=0. 624 linter(:)=1. 625 ! linter(:)=1. 626 b1=2. 627 ! Initialisation des variables entieres 628 lmix(:)=1 629 lmix_bis(:)=2 630 wmaxa(:)=0. 631 lalim(:)=1 632 633 ! print*,'THERMCELL PLUME ARNAUD DEDANS' 634 635 !------------------------------------------------------------------------- 636 ! On ne considere comme actif que les colonnes dont les deux premieres 637 ! couches sont instables. 638 !------------------------------------------------------------------------- 639 active(:)=ztv(:,1)>ztv(:,2) 640 641 !------------------------------------------------------------------------- 642 ! Definition de l'alimentation a l'origine dans thermcell_init 643 !------------------------------------------------------------------------- 644 do l=1,klev-1 645 do ig=1,ngrid 646 if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then 647 alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & 648 & *sqrt(zlev(ig,l+1)) 649 lalim(:)=l+1 650 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) 651 endif 652 enddo 653 enddo 654 do l=1,klev 655 do ig=1,ngrid 656 if (alim_star_tot(ig) > 1.e-10 ) then 657 alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig) 658 endif 659 enddo 660 enddo 661 alim_star_tot(:)=1. 662 663 664 665 !------------------------------------------------------------------------------ 666 ! Calcul dans la premiere couche 667 ! On decide dans cette version que le thermique n'est actif que si la premiere 668 ! couche est instable. 669 ! Pourrait etre change si on veut que le thermiques puisse se déclencher 670 ! dans une couche l>1 671 !------------------------------------------------------------------------------ 672 do ig=1,ngrid 673 ! Le panache va prendre au debut les caracteristiques de l'air contenu 674 ! dans cette couche. 675 if (active(ig)) then 676 ztla(ig,1)=zthl(ig,1) 677 zqta(ig,1)=po(ig,1) 678 zqla(ig,1)=zl(ig,1) 679 !cr: attention, prise en compte de f*(1)=1 680 f_star(ig,2)=alim_star(ig,1) 681 zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & 682 & *(zlev(ig,2)-zlev(ig,1)) & 683 & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) 684 w_est(ig,2)=zw2(ig,2) 685 endif 686 enddo 687 ! 688 689 !============================================================================== 690 !boucle de calcul de la vitesse verticale dans le thermique 691 !============================================================================== 692 do l=2,klev-1 693 !============================================================================== 694 695 696 ! On decide si le thermique est encore actif ou non 697 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test 698 do ig=1,ngrid 699 active(ig)=active(ig) & 700 & .and. zw2(ig,l)>1.e-10 & 701 & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 702 enddo 703 704 705 706 !--------------------------------------------------------------------------- 707 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l 708 ! sans tenir compte du detrainement et de l'entrainement dans cette 709 ! couche 710 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer 711 ! avant) a l'alimentation pour avoir un calcul plus propre 712 !--------------------------------------------------------------------------- 713 714 call thermcell_condens(ngrid,active, & 715 & zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l)) 716 717 718 do ig=1,ngrid 719 ! print*,'active',active(ig),ig,l 720 if(active(ig)) then 721 ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l) 722 zta_est(ig,l)=ztva_est(ig,l) 723 ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) 724 ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & 725 & -zqla_est(ig,l))-zqla_est(ig,l)) 726 727 !------------------------------------------------ 728 !AJAM:nouveau calcul de w² 729 !------------------------------------------------ 730 zdz=zlev(ig,l+1)-zlev(ig,l) 731 zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) 732 733 zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) 734 zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon) 735 w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) 736 737 738 if (w_est(ig,l+1).lt.0.) then 739 w_est(ig,l+1)=zw2(ig,l) 740 endif 741 endif 742 enddo 743 744 745 !------------------------------------------------- 746 !calcul des taux d'entrainement et de detrainement 747 !------------------------------------------------- 748 749 do ig=1,ngrid 750 if (active(ig)) then 751 752 zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) 753 zw2m=w_est(ig,l+1) 754 zdz=zlev(ig,l+1)-zlev(ig,l) 755 zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) 756 ! zbuoybis=zbuoy(ig,l)+RG*0.1/300. 757 zbuoybis=zbuoy(ig,l) 758 zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) 759 zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l) 760 761 762 entr_star(ig,l)=f_star(ig,l)*zdz* zbetalpha*MAX(0., & 763 & afact*zbuoybis/zw2m - fact_epsilon ) 764 765 766 detr_star(ig,l)=f_star(ig,l)*zdz & 767 & *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m & 768 & + 0.012*(zdqt(ig,l)/zw2m)**0.5 ) 769 770 ! En dessous de lalim, on prend le max de alim_star et entr_star pour 771 ! alim_star et 0 sinon 772 if (l.lt.lalim(ig)) then 773 alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) 774 entr_star(ig,l)=0. 775 endif 776 777 ! Calcul du flux montant normalise 778 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & 779 & -detr_star(ig,l) 780 781 endif 782 enddo 783 784 785 !---------------------------------------------------------------------------- 786 !calcul de la vitesse verticale en melangeant Tl et qt du thermique 787 !--------------------------------------------------------------------------- 788 activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 789 do ig=1,ngrid 790 if (activetmp(ig)) then 791 Zsat=.false. 792 ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & 793 & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & 794 & /(f_star(ig,l+1)+detr_star(ig,l)) 795 zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & 796 & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & 797 & /(f_star(ig,l+1)+detr_star(ig,l)) 798 799 endif 800 enddo 801 802 call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l)) 803 804 805 do ig=1,ngrid 806 if (activetmp(ig)) then 807 if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l 808 ! on ecrit de maniere conservative (sat ou non) 809 ! T = Tl +Lv/Cp ql 810 ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) 811 ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l) 812 !on rajoute le calcul de zha pour diagnostiques (temp potentielle) 813 zha(ig,l) = ztva(ig,l) 814 ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & 815 & -zqla(ig,l))-zqla(ig,l)) 816 817 !on ecrit zqsat 818 zqsatth(ig,l)=qsatbef 819 820 zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) 821 zdz=zlev(ig,l+1)-zlev(ig,l) 822 zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) 823 824 zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) 825 zdw2=afact*zbuoy(ig,l)/(fact_epsilon) 826 zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) 827 endif 828 enddo 829 830 if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l 831 ! 832 !--------------------------------------------------------------------------- 833 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 834 !--------------------------------------------------------------------------- 835 836 do ig=1,ngrid 837 if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then 838 ! stop'On tombe sur le cas particulier de thermcell_dry' 839 print*,'On tombe sur le cas particulier de thermcell_plume' 840 zw2(ig,l+1)=0. 841 linter(ig)=l+1 842 endif 843 844 if (zw2(ig,l+1).lt.0.) then 845 linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & 846 & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) 847 zw2(ig,l+1)=0. 848 endif 849 850 wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 851 852 if (wa_moy(ig,l+1).gt.wmaxa(ig)) then 853 ! lmix est le niveau de la couche ou w (wa_moy) est maximum 854 !on rajoute le calcul de lmix_bis 855 if (zqla(ig,l).lt.1.e-10) then 856 lmix_bis(ig)=l+1 857 endif 858 lmix(ig)=l+1 859 wmaxa(ig)=wa_moy(ig,l+1) 860 endif 861 enddo 862 863 !========================================================================= 864 ! FIN DE LA BOUCLE VERTICALE 865 enddo 866 !========================================================================= 867 868 ! print*,'THERMCELL PLUME ARNAUD DEDANS 7' 869 !on recalcule alim_star_tot 870 do ig=1,ngrid 871 alim_star_tot(ig)=0. 872 enddo 873 do ig=1,ngrid 874 do l=1,lalim(ig)-1 875 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) 876 enddo 877 enddo 878 879 880 if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l 881 882 #undef wrgrads_thermcell 883 #ifdef wrgrads_thermcell 884 call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ') 885 call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ') 886 call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ') 887 call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ') 888 call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ') 889 call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ') 890 call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ') 891 #endif 892 893 894 ! print*,'THERMCELL PLUME ARNAUD DEDANS 8' 895 return 896 end 897 -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/yamada4.F
r1299 r1311 1 1 ! 2 ! $ Id$2 ! $Header$ 3 3 ! 4 4 SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp … … 38 38 c iflag_pbl=7 : MY 2.0.Fournier 39 39 c iflag_pbl=8 : MY 2.5 40 c iflag_pbl =9 : un test ?40 c iflag_pbl>=9 : MY 2.5 avec diffusion verticale 41 41 42 42 c....................................................................... … … 66 66 67 67 integer nlay,nlev 68 cym PARAMETER (nlay=klev)69 cym PARAMETER (nlev=klev+1)70 68 71 69 logical first … … 98 96 real fl,zzz,zl0,zq2,zn2 99 97 100 cym real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev) 101 cym s ,lyam(klon,klev),knyam(klon,klev) 102 cym s ,w2yam(klon,klev),t2yam(klon,klev) 103 real,allocatable,save,dimension(:,:) :: rino,smyam,styam,lyam, 104 s knyam,w2yam,t2yam 105 cym common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam 106 c$OMP THREADPRIVATE(rino,smyam,styam,lyam,knyam,w2yam,t2yam) 98 real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev) 99 s ,lyam(klon,klev),knyam(klon,klev) 100 s ,w2yam(klon,klev),t2yam(klon,klev) 107 101 logical,save :: firstcall=.true. 108 109 character (len=20) :: modname='yamada4'110 character (len=80) :: abort_message111 112 102 c$OMP THREADPRIVATE(firstcall) 113 103 frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) … … 123 113 124 114 if (firstcall) then 125 allocate(rino(klon,klev+1),smyam(klon,klev),styam(klon,klev))126 allocate(lyam(klon,klev),knyam(klon,klev))127 allocate(w2yam(klon,klev),t2yam(klon,klev))128 115 allocate(l0(klon)) 129 116 firstcall=.false. … … 131 118 132 119 133 if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then 134 abort_message = 'probleme de coherence dans appel a MY' 135 CALL abort_gcm (modname,abort_message,1) 120 if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then 121 stop'probleme de coherence dans appel a MY' 136 122 endif 137 123 … … 417 403 enddo 418 404 419 c if (iflag_pbl.ge.7..and.0.eq.1) then 420 c q2(:,1)=q2(:,2) 421 c call vdif_q2(dt,g,rconst,plev,temp,kq,q2) 422 c endif 405 if (iflag_pbl.ge.9) then 406 ! print*,'YAMADA VDIF' 407 q2(:,1)=q2(:,2) 408 call vdif_q2(dt,g,rconst,ngrid,plev,temp,kq,q2) 409 endif 423 410 424 411 c Traitement des cas noctrunes avec l'introduction d'une longueur … … 497 484 return 498 485 end 486 SUBROUTINE vdif_q2(timestep,gravity,rconst,ngrid,plev,temp, 487 & kmy,q2) 488 use dimphy 489 IMPLICIT NONE 490 c....................................................................... 491 #include "dimensions.h" 492 cccc#include "dimphy.h" 493 c....................................................................... 494 c 495 c dt : pas de temps 496 497 real plev(klon,klev+1) 498 real temp(klon,klev) 499 real timestep 500 real gravity,rconst 501 real kstar(klon,klev+1),zz 502 real kmy(klon,klev+1) 503 real q2(klon,klev+1) 504 real deltap(klon,klev+1) 505 real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1) 506 integer ngrid 507 508 integer i,k 509 510 ! print*,'RD=',rconst 511 do k=1,klev 512 do i=1,ngrid 513 c test 514 ! print*,'i,k',i,k 515 ! print*,'temp(i,k)=',temp(i,k) 516 ! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1) 517 zz=(plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k)) 518 kstar(i,k)=0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz 519 s /(plev(i,k)-plev(i,k+1))*timestep 520 enddo 521 enddo 522 523 do k=2,klev 524 do i=1,ngrid 525 deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1)) 526 enddo 527 enddo 528 do i=1,ngrid 529 deltap(i,1)=0.5*(plev(i,1)-plev(i,2)) 530 deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1)) 531 denom(i,klev+1)=deltap(i,klev+1)+kstar(i,klev) 532 alpha(i,klev+1)=deltap(i,klev+1)*q2(i,klev+1)/denom(i,klev+1) 533 beta(i,klev+1)=kstar(i,klev)/denom(i,klev+1) 534 enddo 535 536 do k=klev,2,-1 537 do i=1,ngrid 538 denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))* 539 s kstar(i,k)+kstar(i,k-1) 540 c correction d'un bug 10 01 2001 541 alpha(i,k)=(q2(i,k)*deltap(i,k) 542 s +kstar(i,k)*alpha(i,k+1))/denom(i,k) 543 beta(i,k)=kstar(i,k-1)/denom(i,k) 544 enddo 545 enddo 546 547 c Si on recalcule q2(1) 548 if(1.eq.0) then 549 do i=1,ngrid 550 denom(i,1)=deltap(i,1)+(1-beta(i,2))*kstar(i,1) 551 q2(i,1)=(q2(i,1)*deltap(i,1) 552 s +kstar(i,1)*alpha(i,2))/denom(i,1) 553 enddo 554 endif 555 c sinon, on peut sauter cette boucle... 556 557 do k=2,klev+1 558 do i=1,ngrid 559 q2(i,k)=alpha(i,k)+beta(i,k)*q2(i,k-1) 560 enddo 561 enddo 562 563 return 564 end 565 SUBROUTINE vdif_q2e(timestep,gravity,rconst,ngrid, 566 & plev,temp,kmy,q2) 567 use dimphy 568 IMPLICIT NONE 569 c....................................................................... 570 #include "dimensions.h" 571 cccc#include "dimphy.h" 572 c....................................................................... 573 c 574 c dt : pas de temps 575 576 real plev(klon,klev+1) 577 real temp(klon,klev) 578 real timestep 579 real gravity,rconst 580 real kstar(klon,klev+1),zz 581 real kmy(klon,klev+1) 582 real q2(klon,klev+1) 583 real deltap(klon,klev+1) 584 real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1) 585 integer ngrid 586 587 integer i,k 588 589 do k=1,klev 590 do i=1,ngrid 591 zz=(plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k)) 592 kstar(i,k)=0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz 593 s /(plev(i,k)-plev(i,k+1))*timestep 594 enddo 595 enddo 596 597 do k=2,klev 598 do i=1,ngrid 599 deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1)) 600 enddo 601 enddo 602 do i=1,ngrid 603 deltap(i,1)=0.5*(plev(i,1)-plev(i,2)) 604 deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1)) 605 enddo 606 607 do k=klev,2,-1 608 do i=1,ngrid 609 q2(i,k)=q2(i,k)+ 610 s ( kstar(i,k)*(q2(i,k+1)-q2(i,k)) 611 s -kstar(i,k-1)*(q2(i,k)-q2(i,k-1)) ) 612 s /deltap(i,k) 613 enddo 614 enddo 615 616 do i=1,ngrid 617 q2(i,1)=q2(i,1)+ 618 s ( kstar(i,1)*(q2(i,2)-q2(i,1)) 619 s ) 620 s /deltap(i,1) 621 q2(i,klev+1)=q2(i,klev+1)+ 622 s ( 623 s -kstar(i,klev)*(q2(i,klev+1)-q2(i,klev)) ) 624 s /deltap(i,klev+1) 625 enddo 626 627 return 628 end
Note: See TracChangeset
for help on using the changeset viewer.
