Changeset 5852 for LMDZ6


Ignore:
Timestamp:
Oct 10, 2025, 3:57:59 PM (2 months ago)
Author:
fhourdin
Message:

Separation thermcell_plume 5B/6A

Location:
LMDZ6/trunk/libf/phylmd
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume_6A.f90

    r5620 r5852  
    684684 RETURN
    685685     END SUBROUTINE thermcell_plume_6A
    686 
    687 
    688 
    689 
    690 
    691 
    692 
    693 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    694 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    695 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    696 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    697 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    698 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    699  SUBROUTINE thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    700 &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    701 &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
    702 &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    703 &           ,lev_out,lunout1,igout)
    704 !&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
    705 
    706 !--------------------------------------------------------------------------
    707 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    708 ! Version conforme a l'article de Rio et al. 2010.
    709 ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
    710 !--------------------------------------------------------------------------
    711 
    712       USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
    713        USE lmdz_thermcell_qsat, ONLY : thermcell_qsat
    714       IMPLICIT NONE
    715 
    716       INTEGER itap
    717       INTEGER lunout1,igout
    718       INTEGER ngrid,nlay
    719       REAL ptimestep
    720       REAL ztv(ngrid,nlay)
    721       REAL zthl(ngrid,nlay)
    722       REAL, INTENT(IN) :: po(ngrid,nlay)
    723       REAL zl(ngrid,nlay)
    724       REAL rhobarz(ngrid,nlay)
    725       REAL zlev(ngrid,nlay+1)
    726       REAL pplev(ngrid,nlay+1)
    727       REAL pphi(ngrid,nlay)
    728       REAL zpspsk(ngrid,nlay)
    729       REAL alim_star(ngrid,nlay)
    730       REAL f0(ngrid)
    731       INTEGER lalim(ngrid)
    732       integer lev_out                           ! niveau pour les print
    733       integer nbpb
    734    
    735       real alim_star_tot(ngrid)
    736 
    737       REAL ztva(ngrid,nlay)
    738       REAL ztla(ngrid,nlay)
    739       REAL zqla(ngrid,nlay)
    740       REAL zqta(ngrid,nlay)
    741       REAL zha(ngrid,nlay)
    742 
    743       REAL detr_star(ngrid,nlay)
    744       REAL coefc
    745       REAL entr_star(ngrid,nlay)
    746       REAL detr(ngrid,nlay)
    747       REAL entr(ngrid,nlay)
    748 
    749       REAL csc(ngrid,nlay)
    750 
    751       REAL zw2(ngrid,nlay+1)
    752       REAL w_est(ngrid,nlay+1)
    753       REAL f_star(ngrid,nlay+1)
    754       REAL wa_moy(ngrid,nlay+1)
    755 
    756       REAL ztva_est(ngrid,nlay)
    757       REAL zqla_est(ngrid,nlay)
    758       REAL zqsatth(ngrid,nlay)
    759       REAL zta_est(ngrid,nlay)
    760       REAL zbuoyjam(ngrid,nlay)
    761       REAL ztemp(ngrid),zqsat(ngrid)
    762       REAL zdw2
    763       REAL zw2modif
    764       REAL zw2fact
    765       REAL zeps(ngrid,nlay)
    766 
    767       REAL linter(ngrid)
    768       INTEGER lmix(ngrid)
    769       INTEGER lmix_bis(ngrid)
    770       REAL    wmaxa(ngrid)
    771 
    772       INTEGER ig,l,k
    773 
    774       real zdz,zbuoy(ngrid,nlay),zalpha,gamma(ngrid,nlay),zdqt(ngrid,nlay),zw2m
    775       real zbuoybis
    776       real zcor,zdelta,zcvm5,qlbef,zdz2
    777       real betalpha,zbetalpha
    778       real eps, afact
    779       logical Zsat
    780       LOGICAL active(ngrid),activetmp(ngrid)
    781       REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
    782       REAL c2(ngrid,nlay)
    783       Zsat=.false.
    784 ! Initialisation
    785 
    786       fact_epsilon=0.002
    787       betalpha=0.9
    788       afact=2./3.           
    789 
    790       zbetalpha=betalpha/(1.+betalpha)
    791 
    792 
    793 ! Initialisations des variables reeles
    794 if (1==1) then
    795       ztva(:,:)=ztv(:,:)
    796       ztva_est(:,:)=ztva(:,:)
    797       ztla(:,:)=zthl(:,:)
    798       zqta(:,:)=po(:,:)
    799       zha(:,:) = ztva(:,:)
    800 else
    801       ztva(:,:)=0.
    802       ztva_est(:,:)=0.
    803       ztla(:,:)=0.
    804       zqta(:,:)=0.
    805       zha(:,:) =0.
    806 endif
    807 
    808       zqla_est(:,:)=0.
    809       zqsatth(:,:)=0.
    810       zqla(:,:)=0.
    811       detr_star(:,:)=0.
    812       entr_star(:,:)=0.
    813       alim_star(:,:)=0.
    814       alim_star_tot(:)=0.
    815       csc(:,:)=0.
    816       detr(:,:)=0.
    817       entr(:,:)=0.
    818       zw2(:,:)=0.
    819       zbuoy(:,:)=0.
    820       zbuoyjam(:,:)=0.
    821       gamma(:,:)=0.
    822       zeps(:,:)=0.
    823       w_est(:,:)=0.
    824       f_star(:,:)=0.
    825       wa_moy(:,:)=0.
    826       linter(:)=1.
    827 !     linter(:)=1.
    828 ! Initialisation des variables entieres
    829       lmix(:)=1
    830       lmix_bis(:)=2
    831       wmaxa(:)=0.
    832       lalim(:)=1
    833 
    834 
    835 !-------------------------------------------------------------------------
    836 ! On ne considere comme actif que les colonnes dont les deux premieres
    837 ! couches sont instables.
    838 !-------------------------------------------------------------------------
    839       active(:)=ztv(:,1)>ztv(:,2)
    840 
    841 !-------------------------------------------------------------------------
    842 ! Definition de l'alimentation
    843 !-------------------------------------------------------------------------
    844       do l=1,nlay-1
    845          do ig=1,ngrid
    846             if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
    847                alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
    848      &                       *sqrt(zlev(ig,l+1))
    849                lalim(ig)=l+1
    850                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    851             endif
    852          enddo
    853       enddo
    854       do l=1,nlay
    855          do ig=1,ngrid
    856             if (alim_star_tot(ig) > 1.e-10 ) then
    857                alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
    858             endif
    859          enddo
    860       enddo
    861       alim_star_tot(:)=1.
    862 
    863 
    864 
    865 !------------------------------------------------------------------------------
    866 ! Calcul dans la premiere couche
    867 ! On decide dans cette version que le thermique n'est actif que si la premiere
    868 ! couche est instable.
    869 ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
    870 ! dans une couche l>1
    871 !------------------------------------------------------------------------------
    872 do ig=1,ngrid
    873 ! Le panache va prendre au debut les caracteristiques de l'air contenu
    874 ! dans cette couche.
    875     if (active(ig)) then
    876     ztla(ig,1)=zthl(ig,1)
    877     zqta(ig,1)=po(ig,1)
    878     zqla(ig,1)=zl(ig,1)
    879 !cr: attention, prise en compte de f*(1)=1
    880     f_star(ig,2)=alim_star(ig,1)
    881     zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
    882 &                     *(zlev(ig,2)-zlev(ig,1))  &
    883 &                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
    884     w_est(ig,2)=zw2(ig,2)
    885     endif
    886 enddo
    887 !
    888 
    889 !==============================================================================
    890 !boucle de calcul de la vitesse verticale dans le thermique
    891 !==============================================================================
    892 do l=2,nlay-1
    893 !==============================================================================
    894 
    895 
    896 ! On decide si le thermique est encore actif ou non
    897 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
    898     do ig=1,ngrid
    899        active(ig)=active(ig) &
    900 &                 .and. zw2(ig,l)>1.e-10 &
    901 &                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
    902     enddo
    903 
    904 
    905 
    906 !---------------------------------------------------------------------------
    907 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
    908 ! sans tenir compte du detrainement et de l'entrainement dans cette
    909 ! couche
    910 ! C'est a dire qu'on suppose
    911 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
    912 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    913 ! avant) a l'alimentation pour avoir un calcul plus propre
    914 !---------------------------------------------------------------------------
    915 
    916    ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
    917    call thermcell_qsat(ngrid, 1, active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
    918 
    919     do ig=1,ngrid
    920 !       print*,'active',active(ig),ig,l
    921         if(active(ig)) then
    922         zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
    923         ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    924         zta_est(ig,l)=ztva_est(ig,l)
    925         ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
    926         ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    927      &      -zqla_est(ig,l))-zqla_est(ig,l))
    928 
    929 !------------------------------------------------
    930 !AJAM:nouveau calcul de w? 
    931 !------------------------------------------------
    932               zdz=zlev(ig,l+1)-zlev(ig,l)
    933               zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    934 
    935               zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    936               zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
    937               w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
    938  
    939 
    940              if (w_est(ig,l+1).lt.0.) then
    941                 w_est(ig,l+1)=zw2(ig,l)
    942              endif
    943        endif
    944     enddo
    945 
    946 
    947 !-------------------------------------------------
    948 !calcul des taux d'entrainement et de detrainement
    949 !-------------------------------------------------
    950 
    951      do ig=1,ngrid
    952         if (active(ig)) then
    953 
    954           zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
    955           zw2m=w_est(ig,l+1)
    956           zdz=zlev(ig,l+1)-zlev(ig,l)
    957           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    958 !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
    959           zbuoybis=zbuoy(ig,l)
    960           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    961           zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
    962 
    963          
    964           entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
    965     &     afact*zbuoybis/zw2m - fact_epsilon )
    966 
    967 
    968           detr_star(ig,l)=f_star(ig,l)*zdz                        &
    969     &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
    970     &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
    971          
    972 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
    973 ! alim_star et 0 sinon
    974         if (l.lt.lalim(ig)) then
    975           alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
    976           entr_star(ig,l)=0.
    977         endif
    978 
    979 ! Calcul du flux montant normalise
    980       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    981      &              -detr_star(ig,l)
    982 
    983       endif
    984    enddo
    985 
    986 
    987 !----------------------------------------------------------------------------
    988 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
    989 !---------------------------------------------------------------------------
    990    activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
    991    do ig=1,ngrid
    992        if (activetmp(ig)) then
    993            Zsat=.false.
    994            ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
    995      &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
    996      &            /(f_star(ig,l+1)+detr_star(ig,l))
    997            zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
    998      &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
    999      &            /(f_star(ig,l+1)+detr_star(ig,l))
    1000 
    1001         endif
    1002     enddo
    1003 
    1004    ztemp(:)=zpspsk(:,l)*ztla(:,l)
    1005    call thermcell_qsat(ngrid, 1, activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    1006 
    1007    do ig=1,ngrid
    1008       if (activetmp(ig)) then
    1009 ! on ecrit de maniere conservative (sat ou non)
    1010 !          T = Tl +Lv/Cp ql
    1011            zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
    1012            ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
    1013            ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
    1014 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
    1015            zha(ig,l) = ztva(ig,l)
    1016            ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
    1017      &              -zqla(ig,l))-zqla(ig,l))
    1018            zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    1019            zdz=zlev(ig,l+1)-zlev(ig,l)
    1020            zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
    1021 
    1022             zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    1023             zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
    1024             zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
    1025       endif
    1026    enddo
    1027 
    1028         if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
    1029 !
    1030 !---------------------------------------------------------------------------
    1031 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    1032 !---------------------------------------------------------------------------
    1033 
    1034    nbpb=0
    1035    do ig=1,ngrid
    1036             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    1037 !               stop'On tombe sur le cas particulier de thermcell_dry'
    1038 !               print*,'On tombe sur le cas particulier de thermcell_plume'
    1039                 nbpb=nbpb+1
    1040                 zw2(ig,l+1)=0.
    1041                 linter(ig)=l+1
    1042             endif
    1043 
    1044         if (zw2(ig,l+1).lt.0.) then
    1045            linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    1046      &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    1047            zw2(ig,l+1)=0.
    1048         elseif (f_star(ig,l+1).lt.0.) then
    1049            linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
    1050      &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
    1051 !           print*,"linter plume", linter(ig)
    1052            zw2(ig,l+1)=0.
    1053         endif
    1054 
    1055            wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    1056 
    1057         if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
    1058 !   lmix est le niveau de la couche ou w (wa_moy) est maximum
    1059 !on rajoute le calcul de lmix_bis
    1060             if (zqla(ig,l).lt.1.e-10) then
    1061                lmix_bis(ig)=l+1
    1062             endif
    1063             lmix(ig)=l+1
    1064             wmaxa(ig)=wa_moy(ig,l+1)
    1065         endif
    1066    enddo
    1067 
    1068    if (nbpb>0) then
    1069    print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
    1070    endif
    1071 
    1072 !=========================================================================
    1073 ! FIN DE LA BOUCLE VERTICALE
    1074       enddo
    1075 !=========================================================================
    1076 
    1077 !on recalcule alim_star_tot
    1078        do ig=1,ngrid
    1079           alim_star_tot(ig)=0.
    1080        enddo
    1081        do ig=1,ngrid
    1082           do l=1,lalim(ig)-1
    1083           alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    1084           enddo
    1085        enddo
    1086        
    1087 
    1088         if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
    1089 
    1090      return
    1091      END SUBROUTINE thermcell_plume_5B
    1092 END MODULE lmdz_thermcell_plume_6A
Note: See TracChangeset for help on using the changeset viewer.